THE  ANALYTIC  EVALUATION  OF 
MOLECULAR  PROPERTIES  AND  GRADIENTS 
FOR  CORRELATED  WAVEFUNCTIONS 


By 

GARY  WAYNE  TRUCKS 


A DISSERTATION  PRESENTED  TO  THE  GRADUATE  SCHOOL 
OF  THE  UNIVERSITY  OF  FLORIDA  IN  PARTIAL  FULFILLMENT 
OF  THE  REQUIREMENTS  FOR  THE  DEGREE  OF 
DOCTOR  OF  PHILOSOPHY 

UNIVERSITY  OF  FLORIDA 

1988 


P DF  F libraries 


ACKNOWLEDGMENTS 


I thank  my  advisor,  Dr.  Rodney  J.  Bartlett,  for  his  direction, 
knowledge  of  many-body  perturbation  theory  and  coupled-cluster  theory, 
and  my  four-year  research  assistantship.  Despite  an  overwhelming 
schedule,  he  always  found  the  time  to  discuss  my  progress.  I am 
grateful  to  Dr.  George  D.  Purvis  III,  Robert  Harrison,  and  Bill  Laidig 
for  teaching  me  a great  deal  about  computer  programming.  My 
collaborators  for  the  many  different  projects  that  I have  been 
involved  in  over  the  past  five  years  have  also  taught  me  a great  deal; 
Josef  Noga,  John  Watts,  John  Stanton,  Carlos  Sosa  and,  of  course,  Alan 
Salter  especially  come  to  mind. 

I want  to  express  my  deepest  appreciation  to  my  wife,  Teresa.  I 
am  confident  that  the  work  I have  done  over  the  past  five  years  would 
not  have  come  to  pass  without  her  constant  support.  My  parents  and 
grandfather  have  shown  a great  deal  of  faith  in  my  career  decisions, 
and  they  have  encouraged  and  sustained  me  throughout  my  education. 

They  share  a large  part  in  all  that  I have  achieved. 

Gainesville  was  a total  learning  experience  for  me.  I learned  a 
great  deal  about  science  and  even  more  about  life.  Most  of  all,  I 
have  learned  that  the  most  important  things  I can  take  away  from  this 
place  are  the  things  that  I brought  with  me. 


iii 


TABLE  OF  CONTENTS 


Page 

ACKNOWLEDGMENTS iii 

ABSTRACT v 

OVERVIEW 1 

INTRODUCTION 4 

MANY-BODY  THEORY 8 

MANY-BODY  FIRST  DERIVATIVES 28 

Derivation  of  General  Expressions 28 

Elimination  of  Singularities 38 

Identification  of  the  Many-Body  Response  Density 55 

IMPLEMENTATION  OF  MANY-BODY  FIRST  DERIVATIVES 72 

APPLICATIONS  OF  MANY-BODY  FIRST  DERIVATIVES 78 

Stationary  Properties 78 

Molecular  Gradients 82 

CONCLUSION 103 

APPENDIX:  SIMPLIFICATION  OF  THE  FOURTH-ORDER 

FIRST  DERIVATIVE 105 

REFERENCES 113 

BIOGRAPHICAL  SKETCH 117 


iv 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

THE  ANALYTIC  EVALUATION  OF 
MOLECULAR  PROPERTIES  AND  GRADIENTS 
FOR  CORRELATED  UAVEFUNCTIONS 

By 

Gary  Wayne  Trucks 
December,  1988 

Chairman:  Rodney  J.  Bartlett 

Major  Department:  Chemistry 

The  theory  and  implementation  of  analytic  first  derivatives  of  the 
Coupled-Cluster  Doubles  (CCD)  and  full  fourth-order  Many-Body 
Perturbation  Theory  (MBPT(4))  energy  is  presented.  The  analytic 
evaluation  of  one-electron  properties  and  molecular  gradients  is  based 
upon  the  construction  of  a response  density.  This  approach  allows  for 
the  efficient  evaluation  of  all  one-electron  properties  at  a fraction 
of  the  cost  of  traditional  numerical  calculations. 

The  method  is  illustrated  for  both  stationary  and  geometric 
perturbations  for  a variety  of  systems.  A pedagogical  comparison  with 
properties  obtained  by  full  Configuration  Interaction  (Cl)  is  made  for 
CH2,  F”  and  HF.  In  addition,  an  extended  basis  set  for  H2O  is 
investigated.  Analytic  gradient  methods  have  been  used  to  predict 
ground  state  structures,  dipole  moments,  and  harmonic  vibrational 
frequencies.  Examples  for  CCD  gradients  include  H2O,  NH^,  HCN  and 
CH^.  Fourth-order  MBPT  has  been  applied  to  H2O  and  NH^.  The  overall 
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comparison  of  these  methods  with  Cl,  CCSD,  CCSDT-1  and  experiment  are 


very  favorable. 
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OVERVIEW 


Theoretical  investigations  of  atomic  and  molecular  systems  mainly 
involve  finding  approximate  solutions  to  the  Schrodinger  equation. 
Quantum  chemists  usually  employ  the  Born-Oppenheimer  approximation, 
which  gives  rise  to  a Schrodinger  equation  which  involves  an 
electronic  Hamiltonian.  The  most  common  procedure  for  approximately 
solving  the  electronic  problem  is  the  Self-Consistent  Field  (SCF)  or 
Hartree-Fock  (HF)  method.  The  Hartree-Fock  approximation  replaces  the 
many-electron  problem  with  a one-electron  problem  by  defining  an 
effective  potential.  The  effective  potential  allows  for  the  potential 
experienced  by  an  electron  to  be  replaced  with  an  average  potential 
due  to  the  interaction  with  all  other  electrons.  The  solution  of  the 
Hartree-Fock  equations  provides  a set  of  molecular  orbitals  which  form 
a Slater  determinant  that  is  the  best  variational  approximation  to  the 
electronic  wavefunction;  that  is,  the  expectation  value  of  the 
electronic  Hamiltonian  over  the  SCF  wavefunction  provides  the  lowest 
energy  possible  by  a single  determinant  description.  In  practice  the 
Hartree-Fock  equations  are  solved  by  introducing  a finite  set  of 
atomic  basis  functions  (normally  Gaussian  functions  centered  upon  the 
nuclei).  The  molecular  orbitals  are  linear  combinations  of  atomic 
functions.  The  SCF  procedure  is  somewhat  limited  and  is  not  always 
reliable  for  an  adequate  description  of  molecular  structure  and 
properties.  The  reason  is  that  the  repulsion  of  electrons  is  not 
properly  taken  into  account;  a single  determinant  only  introduces 
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exchange  correlation  or  the  motion  of  electrons  with  parallel  spins  is 
correlated  but  the  motion  of  electrons  with  opposite  spins  is  not. 

The  energy  difference  between  the  SCF  and  the  true  nonrelativistic 
energy  is  defined  as  the  correlation  energy.  Correlated  methods 
provide  a means  to  correct  the  single  determinant  and  effective 
potential  approximations  of  SCF  theory.  Indeed,  many  correlated 
methods  build  upon  the  Hartree-Fock.  solution. 

The  Hartree-Fock  determinant  is  only  one  of  the  possible 
determinants  that  may  be  formed  from  the  set  of  molecular  orbitals 
determined  by  the  SCF  procedure.  The  molecular  orbitals  may  be 
labeled  as  occupied  or  virtual  (unoccupied)  orbitals.  Excited 
determinants  may  be  formed  by  mixing  the  SCF  determinant  with  the 
virtual  set  of  orbitals.  For  example,  a singly  excited  determinant 
may  be  formed  by  replacing  an  occupied  orbital  with  a virtual  orbital. 
Higher  excitations  may  be  formed  in  an  analogous  manner.  The  full 
Configuration  Interaction  (Cl)  method  constructs  a matrix 
representation  of  all  possible  determinants.  The  matrix  is 
diagonalized  to  yield  the  energies  of  the  ground  and  excited  states. 
The  Cl  wavefunction  is  the  SCF  determinant  plus  a linear  combination 
of  excited  determinants.  The  full  Cl  wavefunction  is  the  best 
possible  wavefunction  formed  from  a given  basis  set  and  provides  the 
full  correlation  energy. 

In  practice,  the  full  Cl  is  only  available  for  small  systems  due 
to  the  overhead  involved.  If  the  basis  set  consists  of  N basis 
functions,  there  are  approximately  N!  excited  determinants. 

Approximate  Cl  procedures  exist  in  which  the  number  of  excited 
determinants  is  limited  to  a given  excitation  level.  The  most  widely 
used  Cl  approximation  is  the  singles  and  doubles  Cl  (SDCI)  which 
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contains  only  singly-  and  doubly-excited  determinants.  The 
corresponding  vavefunction  is  the  best  variationally  optimized 
wavefunction  in  the  truncated  space. 

The  SCF  wavefunction  may  be  improved  by  allowing  variation  of  the 
orbitals  as  well  as  the  coefficients  of  the  linearly  expanded  set  of 
determinants.  Such  a procedure  exists  and  is  termed  the  Multi- 
Configurational  SCF  (MCSCF)  method. 

Non-variational  methods  may  also  be  used  to  recover  the 
correlation  energy.  Many-body  methods  such  as  Many-Body  Perturbation 
Theory  (MBPT)  and  Coupled-Cluster  (CC)  theory  are  two  such  methods. 
Such  methods  have  proven  to  be  very  accurate  approximations  to  the 
full  Cl.  In  addition,  they  possess  the  size-extensive  property  of  the 
full  Cl  model.  In  order  for  a model  to  be  size-extensive,  the  energy 
must  scale  properly  with  the  size  of  the  molecular  system.  For 
example,  consider  a system  is  composed  of  N non-interacting  identical 
components;  the  full  Cl  energy  will  be  N times  the  energy  of  a single 
component.  This  property  is  not  preserved  in  truncated  Cl  or  MCSCF 
models.  Size-extensivity  is  essential  in  applications  involving  large 
systems,  but  it  has  been  found  to  be  important  even  in  small  systems 
[1].  However,  the  loss  of  the  variational  property  is  not  as 
significant,  since  most  properties  of  interest  are  energy  differences 
which  do  not  retain  the  variational  property  even  if  the  energy  is 
determined  variationally  in  each  case. 


INTRODUCTION 


Energy  derivatives  are  of  key  importance  in  theoretical  studies  of 
atomic  and  molecular  systems.  The  first  derivatives  of  the  energy  with 
respect  to  nuclear  displacement,  termed  molecular  gradients,  may  be 
used  in  conjunction  with  second  derivatives  to  determine  minima  and 
saddle  points  on  the  energy  surface  of  a molecular  system.  Second 
derivatives  give  rise  to  harmonic  frequencies  and  intensities.  In 
addition,  energy  derivatives  may  be  used  to  determine  a variety  of 
molecular  properties,  such  as  the  electric  moments  or  polarizabilities 
given  by  the  change  of  the  energy  in  the  presence  of  an  external 
field,  spin  properties  expressed  in  terms  of  spin  densities  from  which 
hyperfine  coupling  constants  may  be  accessed,  diamagnetic 
susceptibilities  and  chemical  shifts  associated  with  nuclear  magnetic 
resonance  which  are  calculated  as  energy  derivatives  with  respect  to 
external  or  nuclear  magnetic  fields. 

Analytic  derivatives  have  become  commonplace  in  quantum  chemistry. 
The  earlier  work  centered  around  derivatives  of  Hartree-Fock  theory 
[2].  As  the  need  for  more  accurate  methods  became  apparent, 
derivatives  of  correlated  wavefunctions  including  Configuration 
Interaction  (Cl)  [3],  Multiconfiguration  Self-consistent  Field  (MCSCF) 
[4],  and  Second-order  Many  Body  Perturbation  Theory  (MBPT)  [5]  became 
available. 

Derivatives  of  variational  methods  were  developed  prior  to  those 
of  non-variational  methods  due  to  simplifications  that  occur  in  the 
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derivative  expressions.  This  fact  is  illustrated  by  writing  the 
energy  derivative  in  the  form: 

3E(c(X),C(X),I(X))  9E  9c  3E  3C  3E  31 

— = — — + — — + — — 

3X  9c  3X  9C  3X  91  3X 

where  c represents  the  molecular  orbital  coefficients,  C the 
coefficients  in  the  determinantal  expansion  of  the  wavefunction,  and  I 
the  one-  and  two-electron  integrals  involving  the  atomic  orbital  basis 
set  functions.  This  expression  is  simplified  for  any  model  in  which 
the  energy  is  optimum  with  respect  to  one  or  more  of  these  parameters. 
In  the  case  of  the  SCF  model,  the  wavefunction  is  a single 
determinant,  and  it  is  not  necessary  to  evaluate  the  derivative  of  the 
molecular  orbital  coefficients  with  respect  to  the  perturbation 
parameter  (3c/3X)  since  the  energy  is  optimal  with  respect  to  changes 
in  the  molecular  orbital  coefficients  (3E/9c  = 0);  therefore,  only  the 
third  term  need  be  evaluated.  For  the  case  of  the  MCSCF  model,  the 
molecular  orbital  coefficients  and  the  expansion  coefficients  are 
optimized.  Both  3E/3C  and  3E/3c  are  zero,  and  again  only  the  third 
term  remains.  In  Cl  methods,  the  energy  is  optimal  with  respect  to 
the  determinantal  coefficients  (3E/3C  = 0)  only.  This  leads  to  a more 
complicated  expression  for  Cl  energy  derivatives;  nevertheless.  Cl 
energy  derivatives  are  now  routinely  available. 

Although  many-body  perturbation  theory  [6-11]  and  infinite-order 
coupled-cluster  [6,7,8,12-19]  models  have  established  themselves  as  an 
effective  means  to  include  the  effects  of  electron  correlation 
[20,21],  there  are  complications  which  arise  in  the  energy  derivatives 
of  these  methods  not  present  in  variational  approaches.  Due  to  the 
size  extensive  property,  MBPT/CC  methods  are  not  variational;  the 
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energy  is  not  optimal  with  respect  to  either  the  molecular  orbital 
coefficients  or  the  expansion  coefficients.  Therefore,  it  would 
appear  that  the  first-order  change  in  the  wavefunction  must  be 
evaluated.  A direct  treatment  would  require  separate  solutions  of  the 
MBPT/CC  equations  for  each  degree  of  freedom.  The  method  was  applied 
in  this  manner  for  a few  molecular  systems;  however,  the  computational 
expense  prohibited  general  purpose  applications.  However,  the 
equations  that  describe  the  change  in  the  SCF  coefficients  and  the  CC 
expansion  coefficients  are  linear.  In  both  cases  a single 
perturbation-independent  linear  equation  may  be  formulated,  thus 
greatly  simplifying  the  CC/MBPT  gradient. 

This  work  features  an  efficient  formulation  and  implementation  of 
analytical  MBPT/CC  first  derivatives.  First  derivatives  are  necessary 
for  the  expedient  evaluation  of  molecular  structures  in  all  levels  of 
theory;  however,  the  need  for  first  derivatives  in  MBPT/CC  methods  is 
magnified  by  the  fact  that  no  other  efficient  procedure  exists  for  the 
calculation  of  one-electron  properties.  This  work  is  centered  around 
the  identification  and  implementation  of  a perturbation-independent 
response  density  matrix.  The  response  density  allows  one-electron 
properties  to  be  computed  as  dot-products  of  the  response  density  with 
one-electron  property  integrals.  Gradients  are  formulated  in  the  same 
manner  with  the  addition  of  contributions  arising  from  the  derivatives 
of  the  primitive  basis  functions  and  derivatives  arising  from  the 
derivative  of  the  SCF  orthogonality  condition. 

In  this  work,  the  response  density  approach  is  implemented  for  the 
MBPT(2)  through  the  SDTQ-MBPT(4)  models  and  the  CCD  model.  For  each 
level  of  theory,  the  response  density  has  been  identified  and 
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implemented  for  the  first  time.  In  addition,  analytical  molecular 
gradients  for  the  SDTQ-MBPT(4)  model  are  presented  and  implemented  for 
the  first  time.  Recent  developments  in  the  theory  of  coupled-cluster 
derivatives  has  led  to  the  efficient  development  and  implementation  of 
CCSD  and  CCSDT-1  gradients  [22,23].  Although  a subset  of  the  CCSD 
equations,  this  work  presents  an  efficient  computational 
implementation  of  CCD  gradients  for  the  first  time.  Analytic  MBPT(3) 
and  SDQ-MBPT(4)  gradients  have  also  been  recently  reported;  however, 
the  implementation  did  not  exploit  the  linear  nature  of  the  SCF 
derivative  equations.  This  work  also  presents  an  improved  formulation 
and  implementation  of  MBPT(3)  and  SDQ-MBPT(4)  gradient  expressions,  as 
well. 


MANY-BODY  THEORY 


The  time-independent  Schrodinger  equation  for  a molecular  system 
is  given  by 


H Y = E T (1) 


where  H is  the  Hamiltonian  operator  for  a system  comprising  electrons 
and  nuclei,  the  wavefunction  Y is  an  eigenfunction  of  H and  the 
eigenvalue  E is  the  energy  of  the  system.  In  atomic  units,  the 
Hamiltonian  is  divided  into  five  parts 


n N n N 

12  1 12 

H.-E  _v__E  V-EEZ/r 

2 i 2 2M  a a ia 

i=l  0=1  “ i=l  0=1 

n N 

+E  1/r  +EZZ/R 

ij  o 3 og 

i>j  <x>3 


where  the  labels  i,j,n  refer  to  electrons  and  o, 3,N  refers  to  nuclei, 

M is  mass  for  the  nuclei,  and  r and  R refer  to  the  distance  between 
o 

particles.  The  first  two  terms  are  the  kinetic  energy  operator  for 
electrons  and  nuclei,  respectively.  The  third  term  is  the  coulomb 
attraction  between  electrons  and  nuclei.  The  repulsion  between 
electrons  and  between  nuclei  are  represented  by  the  fourth  and  fifth 
terms,  respectively. 
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Due  to  the  greater  mass  of  the  nuclei,  they  move  much  more  slowly 
than  electrons;  in  fact,  the  electrons  may  be  thought  of  as  moving  in 
a field  of  fixed  nuclei.  This  is  the  Born-Oppenheimer  approximation 
which  is  central  to  the  study  of  electronic  structure.  This 
approximation  allows  the  kinetic  energy  of  the  nuclei  to  be  neglected; 
the  nuclei  repulsion  is  treated  as  a constant.  The  remaining  terms  of 
the  Hamiltonian  define  the  electronic  problem 


n 

1 

H - - E - 
2 

i=l 


2 

7 

i 


n N n 

-EEZ/r  +El/r 

a ia  ij 

i=l  ot=l  i>j 


The  Schrodinger  equation  becomes 

H_Y_  = E_Y 
e e e e 


The  electronic  problem  only  parametrically  depends  on  the  nuclear 
coordinates  through  the  orbital  basis  functions.  The  total  energy  of 
the  system  is  defined  by  adding  the  constant  nuclear  repulsion  term. 
The  electronic  problem  will  be  the  focus  of  the  remainder  of  this 
work. 

The  Hartree-Fock  procedure  mentioned  earlier  is  an  approximate 
solution  of  the  electronic  problem.  Correlated  methods  provide  more 
accurate  solutions  of  the  electronic  problem.  Correlated  methods 
normally  build  upon  some  reference  function,  i.  In  many  cases  the 
reference  function  is  the  single-determinant  wavefunction  defined  by 
the  Hartree-Fock  solution,  although  others  may  be  chosen.  The  energy 
of  the  reference  function  is  given  as  the  expectation  value 
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E = $ H $ dv 

ref 


or  more  conveniently  expressed 

E = <01H|0>  . 

ref 

It  is  now  common  to  express  operators  in  the  so-called  second 
quantized  or  occupation  number  form.  There  exists  a one-to-one 
correspondence  with  the  ordinary  or  "first-quantized"  form.  The 
Hamiltonian  may  be  expressed 


+ ^ 
H=Ihpq+l/4  E <pq||rs>p  q sr  . 
p,q  pq  P>q>r»s 

By  convention,  the  indices  i,j,k,l,...  represent  the  spin  orbitals 
which  are  occupied  in  the  reference  function  |0>.  The  indices 
a,b,c,d,...  correspond  to  the  virtual  (or  unoccupied)  spin  orbitals. 

The  indices  p,q,r,s,...  are  unrestricted. 

The  symbols  hp^  and  <pql|rs>  in  the  previous  equation  represent 
integrals  wherein  the  one-  and  two-electron  operators  of  the 
Hamiltonian  act  upon  the  orbital  basis  functions: 


' * 

f 1 2 -1  ) 

h = 

pq 

« 

♦p(i) 

7,  - E Z r T 

( Z ot 

<t>q(l) 


* 

' * 

* 

■ -1 

<pq 1 1 rs>  = 

♦p(l) 

♦q(2) 

^12  ^12^ 

4>s(2) 


The  variables  and  are  the  coordinates  of  spin  and  space  of 
electron  number  one,  respectively.  The  operator  P^2 
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permutation  operator  which  Insures  that  the  integral  <pq||rs>  is 
antisymmetric  with  respect  to  the  exchange  of  p with  q or  r with  s. 

The  operators  p"*^  and  q are  creation  and  annihilation  operators, 
respectively.  These  operators  create  excited  determinants  by  acting 
on  the  reference  function.  An  electron  is  either  annihilated  from  an 
occupied  orbital  or  created  in  a previously  unoccupied  or  virtual 
orbital.  This  is  demonstrated  in  the  following  examples: 


{a^} 

|0>  = 

|S> 

Singly-excited 

determinant 

{a'^ib'^j} 

|0>  = 

|D> 

Doubly-excited 

determinant 

(a’^ib’^jc'^k)  |0>  = 

1T> 

Triply-exci ted 

determinant 

In  the  previous  equation,  the  brackets  enclosing  a second- 
quantized  operator  string  {ABCD...}  denote  the  normal  ordered  form  of 
operators.  The  normal  ordered  form  is  defined  as  the  product  which 
has  all  hole  creation  operators  or  particle  annihilation  operators 
(!■*■,  .. . and  a,b,c,...)  to  the  right  of  the  other  operators.  This 

may  be  accomplished  by  commuting  the  operators  such  that  {ABCD...}  = 
{ABDC...}.  This  is  done  to  exploit  the  fact  that  i'^|0>  = 0 and  a|0> 
0 (i.e.  an  orbital  that  is  previously  occupied  cannot  be  doubly 
occupied  and  an  unoccupied  orbital  cannot  be  annihilated).  In  fact, 
these  conditions  illustrate  the  fundamental  anticommutator  relation: 

+ + 

[P,q]^  = p q + qp  = • 

If  non-orthogonal  functions  are  used,  the  overlap  integral  <p|q> 

replaces  S . Furthermore,  a contraction  between  operators  p and  q 

pq 
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may  be  defined  as 

n 

pq  = pq  - {pq}  • 

Application  of  the  anticommutator  relation  allows  a contraction  to 
only  be  nonvanishing  when 

rri  rn  ^ 

i j = and  ab  = . 

At  this  point,  Wick's  Theorem  [25]  may  be  applied.  The  theorem 
states  that  a string  of  second  quantized  operators  is  equal  to  the 
normal  product  of  the  string  plus  all  possible  contractions. 


ABCD...  = {ABCD...}  + {ABCD...}  + (ABCD...)  + {Ako...}  + {AkD...} 

+ all  other  single  contractions 


+ {ABCD...}  + (ABCD...)  + {ABC^..}  + {AkD...} 


+ all  other  double  contractions 


+ all  full  contractions 

Of  course,  all  contractions  within  a normal  ordered  product 
vanish.  There  exists  a Generalized  Wick's  Theorem  [25]  for  operators 
in  normal  ordered  form  that  states  a product  of  a string  of  normal 
ordered  operators  is  equal  to  the  normal  product  of  the  overall  string 
plus  the  sum  of  all  possible  contractions;  contractions  between  pairs 
of  operators  within  an  original  normal  product  are  not  allowed.  A 
corollary  is  that  for  the  expectation  value  of  a product  of  normal 
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product  operators  over  the  reference  (referred  to  as  the  Fermi 
vacuum),  only  the  fully  contracted  products  are  non-zero. 

The  explicit  form  of  the  normal  product  Hamiltonian  in  the  second- 
quantized  or  occupation  number  representation  is  given  by 


h + E <pi I |qi> 

pq 


p»q 


{p'^q}  + - E <pq  I I rs>{p'^q'^sr}  . (2) 

4 

p»q» 


r,s 


Subtraction  of  the  reference  energy  from  the  electronic 
Schrodinger  equation  yields 

(H  - <0|H|0>)  Y = (E  - <0|H10>)  Y 
or 

H.,  Y = AE  Y . 

N 


(3) 


In  the  above  equation,  the  normal  product  form  of  the  Hamiltonian  is 
introduced.  Application  of  Wick's  theorem  results  in  the  following 
normal  ordered  form  of  the  Hamiltonian 


1 

4 


" E fpq{p"^q}  + - <pq|  |rs>{p'^q‘^sr} 


(4) 


p»q 


p.qi 


r.s 


or 


where  f is  the  normal  Fock  operator.  The  representation  of  the  Fock 

pq 

operator  in  the  orbital  basis  is  given  by  the  matrix  elements 


f = h + E <pi 1 I qi>  . 

pq  pq 


i 
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If  the  SCF  wavefunction  is  chosen  as  the  reference  function, 

<0|H|0>  is  the  SCF  energy,  AE  is  defined  as  the  correlation  energy, 
and  the  Hamiltonian  simplifies  to 

Hjj(X)  = E ^ ^ ^ 4 ^ <pq|  |rs>{p‘'q^sr}  , 

i,j  a,b  P,q, 

r,s 

since  the  matrix  elements  f . and  f.  are  zero  for  SCF  orbitals.  For  the 

0.1  Id 

particular  choice  of  canonical  SCF  orbitals,  f.^  = j ^nd  f^j^  = S^^^a’ 

where  e denotes  the  molecular  orbital  eigenvalue  of  orbital  p. 

P 

Many-body  methods  are  used  to  approximate  the  solution  to  the 
Schrodinger  equation  in  this  investigation.  Many-body  perturbation 
theory  (MBPT)  [6-11]  and  coupled  cluster  (CC)  theory  [6,7,8,12-19]  are 
among  the  most  common  many-body  methods.  These  methods  are  size- 
extensive,  meaning  that  the  energy  scales  properly  with  the  size  of 
the  system.  The  methods  were  originally  developed  for  use  in  nuclear 
and  solid  state  physics  where  the  necessity  for  models  that  scaled 
properly  the  the  size  of  the  system  was  crucial.  Size-extensivity 
also  has  important  ramifications  in  chemical  applications  involving 
large  molecules  or  dissociation. 

Although  MBPT  was  developed  prior  to  CC  theory,  there  is  an 
important  relationship  between  the  two  methods.  Coupled-cluster 
theory  is  the  infinite-order  extension  of  perturbation  theory.  In 
fact,  many-body  perturbation  theory  is  generated  by  the  low  order 
iterations  of  coupled-cluster  theory.  Therefore,  CC  theory  will  be 
developed  first,  with  various  orders  of  perturbation  theory  being 
obtained  from  the  CC  equations. 


15 


Coupled-cluster  theory  may  be  formally  written  as  an  exact 
solution  of  the  Schrodinger  equation  by  defining  the  wavefunction  as 
an  exponential 


H Y = E Y 


Y = e |0> 


(5) 


where  T is  defined  by  a sum  of  excitation  operators 


T = + T2  + T3  + . . . + 


The  excitation  operators  are  given  by 


= E t {a'^i} 
i 

i,a 

1 ab 

= - E t {a'^ib'^j} 

4 ij 

i»j 

a,b 


N 


ab. . 


(N!)^  ij 

1 > j > • • • > N 

a^by . • . jN 


{a'^ib'^j  . . . } , 


where  the  t-amplitudes  are  antisymmetrized  coefficients  to  be 
determined.  Substituting  into  the  Schrodinger  equation  yields 


e'^|0>  = AE  e'^|0> 
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T 

Multiplying  on  the  left  by  e“  , the  Schrodlnger  equation  becomes 

e~'^H„  e’^10>  = AE|0>  . (6) 

N 

The  Schrodlnger  equation  may  be  simplified  by  use  of  the  Hausdroff 
commutator  expansion 

e-'^Hjje'^  = + [H^,T]  + 1/2  [[Hj^,T],T]  + 1/3!  [ [ ,T]  ,T]  + 

1/4!  l[[[Hj^,T],T],T],T]  . 

The  expansion  terminates  after  the  four-fold  commutator  due  to  the 
fact  the  is  at  most  a two-particle  operator.  The  products  involved 
may  be  classified  as  either  connected  or  disconnected.  A connected 
part  is  defined  by  a contraction  between  factors  of  the  product  or, 
simply,  there  exists  one  or  more  common  indices  in  the  summation. 
Likewise,  a disconnected  part  has  no  contractions  between  factors, 
i.e.  the  sums  may  be  carried  out  independently.  The  disconnected 
parts  of  the  commutator  cancel,  leaving  only  connected  parts.  As  an 
example,  consider  the  disconnected  part  of  the  commutator  of  [Hj^,T]: 

^^N’'^Jdisconnected  = V " 


+ 

ab  + + ' 

or 

Eh  {p  q} 

E t {a  ib  j} 

1 pq  J 

1 ij 

P,q  i»j 


a,b 


ab  + + ' 

f + 

E t {a  ib  j} 

Eh  {p  q} 

1 ij  J 

1 pq 

\ 

i.j 


a,b 


p.q 
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where  p,q  i,j,a,b.  Furthermore,  if  the  left  factor  of  a product  is 
an  excitation  operator,  there  can  be  no  connected  contribution. 
Consequently,  the  Schrodinger  equation  may  be  written: 

e'^)^|0>  = AE|0>  , (7) 

where  c indicates  the  limitation  to  connected  diagrams. 

By  projecting  Equation  (7)  on  the  left  by  <0|,  we  obtain  the  CC 
energy  expression: 


AE  = <0|(H^  e'^)^10>  . (8) 

The  t-amplitudes  are  determined  by  solving  the  set  of  coupled 
equations  obtained  by  left  projection  of  the  singly-,  doubly-, triply- 
...  excited  determinants  (S,  D,  T,  ...): 

<S|(H^  e^)^|0>  = 0 (9) 

<D|(Hj^  e'^)^|0>  = 0 

<T|(H^  e'^)^|0>  = 0 


In  order  to  evaluate  the  coupled-cluster  energy  the  exponential 
operator  is  expanded  and  fully  contracted  terms  are  identified.  The 


nonvanishing  terms  are 
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2 

AE  = <0|  ( T + T + 1/2  T ) |0>  . 

12  1 c 


1 c 


This  may  be  represented  algebraically  by 


a 


ab  a b b a 


AE  = E f t + 1/4  E <ij 1 tab>  ( t 
i,a  ia  i i,j 


+ t t - t t ) . 

ij  i j 1 j 


a,  b 


Nothing  else  contributes  to  the  energy.  All  other  clusters  may  only 
contribute  indirectly  by  their  contribution  to  t^  and  t2* 

It  should  be  noted  that  no  new  approximation  has  been  introduced 
at  this  point.  In  the  limit  in  which  all  cluster  operators  are 
included,  the  full  Cl  wavefunction  and  energy  are  recovered.  If  the 
expansion  of  T is  truncated,  various  coupled-cluster  approximations 
are  obtained.  The  coupled-cluster  doubles  model  (CCD)  is  defined  if  T 
= T2,  the  CCSD  model  [19]  is  defined  by  T = Tj^  + T2  and  the  CCSDT 
model  is  defined  by  T = T^  + T2  + T^.  Other  approximate  models  may  be 
defined  by  restricting  a cluster  operator  to  linear  contributions. 

The  LCCD  model  is  obtained  if  only  linear  T2  contributions  are 
considered,  the  CCSDT-1  model  [18]  is  defined  by  exp(T)  = exp(Tj^  + T2) 
+ • 

It  is  convenient  make  use  of  diagrammatic  techniques  [27]  to 
identify  all  of  the  fully  contracted  products  in  the  CC  energy  and 
amplitude  equations.  The  diagrams  are  then  interpreted  algebraically 
by  a simple  set  of  rules.  The  operators  f^^  and  (the  normal  product 
Hamiltonian  Hj^)  and  the  coupled-cluster  excitation  operators  Tj^,  T2, 
and  T^  are  represented  diagrammatically  in  Figure  (1)  and  Figure  (2), 
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respectively.  The  diagrams  are  constructed  by  connecting  the 
diagrammatic  representation  of  operators  together  to  form  a given 
excitation  level.  The  excitation  level  of  each  diagram  is  given  by 
the  number  of  pairs  of  open  lines  at  the  top  of  the  diagram.  Zero 
excitations  define  the  CC  energy  and  are  diagrammed  in  Figure  (3). 

The  diagrams  for  the  T^,  T2  and  T^  amplitudes  of  the  CCSDT  model  are 
shown  in  Figures  (4),  (5),  and  (6),  and  have  excitation  levels  of 
one,  two  and  three,  respectively.  The  explicit  algebraic  expressions 
for  the  CCSDT  model  have  been  presented  elsewhere  [18,19]. 

In  order  to  make  the  connection  of  MBPT  and  CC,  it  is  convenient 
to  make  the  choice  of  canonical  orbitals  as  defined  by  an  SCF 
reference  determinant.  For  an  SCF  reference,  only  the  diagonal  parts 
of  the  terms  involving  the  f^^  operator  (i.e.  the  orbital  eigenvalues) 
are  nonzero.  These  terms  may  be  taken  to  the  other  side  of  cluster 
equations,  thus  defining  an  equation  that  may  be  solved  iteratively. 
This  involves  diagrams  (3a, b)  of  figure  (4)  for  the  T^  equation, 
diagrams  (2a, b)  in  figure  (5)  for  the  T2  equation,  and  diagrams  (3a, b) 
of  figure  (6)  for  the  T^  equation.  The  resulting  equations  take  the 
form; 


a 


(E  - C ) t 

i a i 


F (I,T) 
1 


ab 


(e+E-e-e)t  =F  (I,T) 
i j a b ij  2 


(e+e+e-e-e 
i j k a b 


e ) t 


c 


abc 

= F (I,T) 
i jk  3 
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where  Fjj(I,T)  denotes  the  remaining  expression  is  a function  of 
integrals  and  t-amplitudes. 

The  relationship  of  the  CC  and  MBPT  approaches  become  clear  if  one 
initially  focuses  on  the  equations.  If  in  the  right  hand  side  of 
the  equation  is  set  to  zero,  the  first  iteration  of  the  T2  equation 
is  given  by  diagram  (1)  of  figure  (5).  The  justification  for  this 
choice  is  that  the  terms  containing  t2  on  the  right  are  assumed  to  be 
small  compared  to  the  left  hand  side  of  the  equation.  The  resulting 
amplitude  expression  is  obtained; 

ab 

t (1)  = <ab| |ij>  / (e  + £ - e - e ) • 

ij  i j a b 

Inserting  this  value  into  the  the  CC  energy  expression  as  given  in 
Equation  (8)  yields  the  second-order  MBPT  energy.  The  first 
contribution  to  MBPT  is  second  order  due  to  the  fact  that  the  SCF 
energy  is  complete  through  first-order.  Notice  that  the  MBPT(2) 
energy  originated  from  doubly-excited  determinants. 

A second  iteration  may  be  carried  out  by  substituting  the 
previously  obtained  t-amplitude  into  the  linear  terms  of  the  T2 
equation  as  given  by  diagrams  (2c-e)  in  Figure  (5),  and  dividing  by 
the  eigenvalue  term  on  the  left  hand  side  of  the  T2  equation.  The 
third-order  MBPT  energy  is  obtained  by  substitution  of  this  result 
into  equation  (8).  Again,  the  MBPT  energy  is  obtained  from  doubly- 
excited  determinants. 

In  order  to  obtain  the  MBPT(4)  energy,  the  new  t-amplitude  is 
substituted  in  the  linear  terms  and  the  quadratic  terms  of  the  T2 
equation  given  by  diagrams  (2c-e)  and  (5a-d)  of  Figure  (5), 
respectively.  This  defines  only  the  fourth-order  doubles,  D-MBPT(4), 
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and  quadruples,  Q-MBPT(4).  The  fourth-order  quadruples  arise  from  the 

1/2  exponential  wavefunction  and  correspond  to  two 

simultaneous  interactions  of  two  electrons  (a  pair  of  two  electron 

clusters).  The  true  four-body  excitations  occur  first  in  fifth-order 

via  T^,  thus  indicating  that  pairs  of  two-body  clusters  interacting 

are  more  important  [28].  The  remaining  fourth-order  contributions 

come  from  the  initial  contributions  of  T^  and  T^  to  the  T2  equations. 

3.b 

The  initial  iteration  of  T^^  is  obtained  by  substitution  of 
diagrams  (la,b)  of  Figure  (6)  and  dividing  by  the  eigenvalue  term  (e^ 

- c ).  This  results  in  the  fourth-order  singles,  S-MBPT(4),  arising 

Sl 

from  diagrams  (3a, b)  of  Figure  (5).  Likewise,  the  fourth-order 

3 b 

triples,  T-MBPT(4)  result  from  inclusion  of  into  diagrams 

(la,b)  of  Figure  (6),  dividing  by  the  eigenvalue  sum  (e^  + - 

e - e.  - c ),  and  substitution  into  diagrams  (4b, c)  of  Figure  (5). 
a b c 

No  higher  excitation  may  be  formed  in  fourth-order  [29]. 

The  coupled-cluster  equations  have  been  implemented  through  CCSDT 
[19];  however,  CCSD  [6]  and  CCSDT-1  [10]  make  up  the  majority  of 
general  purpose  applications.  The  reason  for  this  is  that  the 
asympotic  dependence  of  the  CCSD  equations  scale  as  where  N is  the 
number  of  basis  functions.  Any  introduction  of  triples  scales  as  N^, 
thus  greatly  increasing  computational  overhead  especially  for  an 
iterative  procedure.  The  current  state  of  the  art  for  perturbation 
theory  is  the  full  SDTQ-MBPT(4)  model  [29].  Again,  basis  set 
dependence  is  due  to  the  inclusion  of  triples;  however,  the 
popularity  of  MBPT  rests  upon  the  fact  that  it  benefits  from  the  size- 
consistent  property  of  coupled-cluster,  yet  it  is  not  solved 
iteratively,  thus  allowing  a relatively  inexpensive  approach  to 
systemmatically  introduce  correlation. 
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The  CC  and  MBPT  methods  have  favorably  compared  to  the  energy 
obtained  by  the  full  Cl  [20,21].  It  should  be  apparent  that  coupled- 
cluster  is  superior  to  truncated  Cl  due  to  the  higher  excitations 
introduced  by  the  exponential  ansatz.  The  coupled-cluster  equations 
are  connected,  thereby  insuring  that  the  resulting  energy  expressions 
are  connected  and  closed  (i.e.  no  net  excitation);  therefore,  the  CC 
energy  is  linked  (i.e.  there  are  no  closed  disconnected  terms).  The 
result  of  a linked  energy  expression  gives  rise  to  the  size-extensive 
property  of  MBPT/CC  models. 


fpq  { P^P}  = 


-X  + 


+ X + 


<pq||rs>  {ptq+sr}  = 


+ 


A 

Jv  ^ 


+ 


* A'A 

A 

A 


Figure  1.  The  Diagrammatic  Representation  of 


t°{Q^i}=  V 

t“^{a^,bVk}=  VVV 


Figure  2.  The  Diagrammatic  Representation  of  Coupled-Cluster 
Excitation  Operators  Tj^,  T2  and  T^. 
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Figure  3.  The  Diagrammatic  Representation  of  the  Coupled-Cluster 
Energy. 
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Figure  4. 


The  Diagrammatic  Representation  of  the  T^^  Amplitude 
Equation  in  the  CCSDT  Model. 
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Figure  5.  The  Diagrammatic  Representation  of  the  T2  Amplitude 
Equation  in  the  CCSDT  Model. 
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Figure  6.  The  Diagrammatic  Representation  of  the  Amplitude 

Equation  in  the  CCSDT  Model. 
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Figure  6 — continued 


MANY-BODY  FIRST  DERIVATIVES 


Derivation  of  General  Expressions 


Consider  the  solution  of  the  Schrodlnger  Equation  where  some 
perturbation  Q has  been  introduced.  The  perturbation  may  be  an  atomic 
displacement  or  an  external  one-electron  potential  energy 
perturbation.  The  strength  of  the  perturbation  is  dependent  upon  some 
parameter  X>  therefore,  the  electronic  Hamiltonian  H(X),  its 
eigenfunction  t(X),  and  its  eigenvalue  E(X)  are  dependent  upon  the 
parameter  X* 

H(X)  Y(X)  = E(X)  T(X)  (10) 

2 -1  -1 
H(X)  - -1/2  £ Vj  - E Z„r^i  * ? % * ' 

i i,a  1 


The  total  energy  of  the  system  in  the  presence  of  the  perturbation 
may,  for  small  X»  be  expanded  in  a power  series 

E(X)  = Eq  + XE^  + 1/2  XE2X  + ••• 


where 


El  = 


3E(X) 

3X 


X=0 


3^E(X) 

3X3X 


X=0  , 


The  properties  evaluated  in  this  manner  are  termed  "response 
properties"  since  the  system  is  allowed  to  respond  to  the 
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perturbation.  If  the  perturbation  is  an  external  electric  field,  for 
example,  is  the  permanent  dipole  moment  of  the  unperturbed  system, 
E2  is  the  polarizability,  and  E^,  E^,  etc.  are  higher 
hyperpolarizabilities.  It  should  be  clear  that  response  properties 
may  be  computed  numerically  or  analytically.  In  this  work,  the 
analytic  route  is  developed,  but  first  it  is  proper  to  consider  the 
numerical  approach. 

First-  and  second-order  response  properties  may  be  evaluated  by 
calculating  the  total  electronic  energy  with  the  perturbation  present 
for  both  plus  and  minus  values  of  X-  This  allows  for  the  cancellation 
of  higher  powers  of  X«  This  method  is  termed  the  "finite-field" 
approach  [30,31].  The  advantage  of  the  approach  is  its  ease  of 
implementation  allowing  for  the  determination  of  properties  with 
existing  computer  codes.  This  method  has  proven  to  be  very  accurate 
provided  the  proper  value  of  the  field  strength  (X)  is  chosen; 
unfortunately,  no  prescription  for  determining  the  proper  field 
strength  is  available,  short  of  trial  and  error. 

The  response  treatment  of  the  coupled-cluster  energy  is  expressed 
as: 


Hj^(X)  e'^^^^|0>  = AE(X)|0>  (H) 

where  the  normal  ordered  form  of  the  Hamiltonian  is  used.  This 
implies  that  the  coupled-cluster  amplitude  equations  were  solved  in 
the  presence  of  the  perturbation.  In  fact,  the  finite  field  approach 
has  been  the  traditional  method  for  computing  many-body  properties 
other  than  the  energy. 
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In  the  previous  section,  the  zeroth-order  equations  for  many-body 
theory  were  developed.  These  are,  of  course,  the  original  equations 
in  the  absence  of  the  perturbation  (X=0). 

<0|  (H^  e'^)j0>  = AE 

and 

<♦1  (Hfj  e'^)jO>  = 0 


T = T(X) 


X=0 


X=0 


AE  = AE(X) 


x=o 


Analogous  quantities  may  be  defined  for  the  perturbed  case 


T 


X 


3T(X) 


9X 


X=0 


N 


X _ 
N " 


'5^ 

. 9X 


9Vf^(X)' 


9X 


x=o 


AE^ 


9AE(X) 


9X 


X=0 


The  treatment  of  the  coupled-cluster  response  equations  has  been 
given  by  many  authors  [22-24,32-39].  A new  approach  termed  the 
"density  approach"  to  gradients  has  been  presented  in  our  earlier  work 
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and  the  theory  has  been  generalized  to  first  and  second  derivatives  of 
coupled-cluster  theory  by  Salter  [40].  In  this  work,  the  "density 
approach"  will  be  developed  and  implemented  for  the  coupled-cluster 
doubles  (CCD)  and  SDTQ-MBPT(A)  models  of  theory. 

The  straightforward  derivative  of  Equation  (10)  with  respect  to  X» 
evaluated  at  X-0,  Is  given  by  the  following,  where  the  fact  that 
excitation  operators  commute  with  one  another  has  been  utilized: 


[(H,  e"),,  tX|  . <hX 


|0> 


Or  more  simply. 


AE^|0>  . 


((H„  e"^)  T^)  + 

" N c c 


<"n 


|o> 


AE^|0> 


The  symbol  |P>  is  defined  to  represent  the  space  of  orthonormal 
determinants  in  which  the  amplitudes  of  T(X)  and  AE(X)  are  determined. 
Furthermore,  define  the  symbol  | $>  to  represent  all  excited 
determinants 


1 = |PXP|  (12) 
|P>  = |0>  + |$>  . 

Projection  on  the  left  by  <P|  onto  the  derivative  equation  yields  an 
expression  for  and  a set  of  linear  equations  which  determine  the 
derivative  t-amplitudes: 


<0|((Hjj  e'^)^T^)^|0>  + <0l(H^  e'^)^|0>  = AE^  (13) 


and 
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<*|((H^  e'‘’)J^)^|0>  + e'^)^|0>  = 0 . (U) 

If  the  unperturbed  coupled-cluster  amplitudes  are  available,  the 
derivative  amplitude  equations  can  be  constructed  and  solved.  The 
derivative  energy  AE^  can  then  be  evaluated  by  Equation  (13).  This 
method  was  used  in  early  reports  of  coupled-cluster  gradients  [35,36]; 
however,  it  is  not  practical  for  general  use. 

The  construction  of  derivative  amplitudes  is  costly  and 
fortunately  can  be  avoided.  If  the  resolvent  1 = |PXP|  may  be 
inserted  into  the  terms  involving  T^  and  the  result  is  left  projected 
by  <P|,  one  obtains 

<0|-T^lPXPl(H^e'^)^|0>  + <0|(Hj5e'’^)^|0>  + <0  | | PXP  |T^|  0>  = AE^ 

and 

<*|-T^lPXP|(Hj^e'^)^|0>  + <*|(Hje'^)^|0>  + <t|  (H^e^)^|PXP|T^|0>  = 0 . 

Several  simplifications  occur  in  the  above  equations  due  to  the  fact 

that  T is  the  solution  of  the  unperturbed  equation  <#|(Hj^e  )^|0>=0  and 

<0l(H„e'^)  |0>*AE.  Also,  note  that  the  term  <0|T^lP>  cannot  be  nonzero 
' N C 

since  T^  is  an  excitation  operator.  The  equations  simplify  to  the 
following 


<0|(hJ  e'^)^|0>  + <0|(Hj^  e^)^|$X$|T^lO>  = AE^ 

and 

<*1-T^|0>AE  + <*|(hJ  e^)^|0>  + e'^)^ | $X$|T^| 0>  = 0 . 
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Y 

If  the  second  equation  is  solved  for  <$|T  |0>; 

<*1T^10>  = - <i\  AE  | $>  <#|(hJ  e'^)^|0>  , 

and  then  substituted  for  <*|T^|0>  in  the  AE^  expression. 

<0|(HjJe'^)^|0>  - <0|(Hj^e'^)^|*>  <$|  (H^e^)^-AE  j $>  <$|(hV^^|0>  = AE^  . 

Now  define  an  antisymmetrized,  perturbation-independent,  de- 
excitation operator  A.  Let  A = + A2  + A^  +...,  where 

i 

A^  «=  £ X { i a}  , 

a 

i,a 

1 ij 

A.,  = - £ X {i'^aj'^b)  , 

^4  ab 

i.j 

a,  b 


1 ij • • • ^ 

A = £ X {iajb. ..},•••  > 

(n! ) ab. . . 

1 , J , • • • 

a , b , . . . 

For  every  T in  T there  is  a corresponding  A in  A.  That  is,  for  the 
n u 

CCD  model  ^ For  the  CCSD  model  A = A^  + A^,  and  for  the  CCSDT 

model  A = A^  + A2  + A^.  The  X-amplitudes  are  required  to  satisfy  the 


condition: 
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<0|A|*>  = - <0|(Hj^  e'^)^|$>  <i\  (Hj^  e^)^-  AE  | $> 

y X 

The  resulting  expression  for  AE'^  is  not  dependent  upon  T . 

<0|(hJ  e^)^|0>  + <0|A|$X$|(hJ  e^)^|0>  = AE^  , 

which  can  be  evaluated  by  solving 

<0lAl*X*|  (Hj^  e’’^)^-  AE  |*>  + <0|(Hj^  e'^)^|$>  = 0 . 

This  method,  referred  to  as  the  Z-vector  method,  was  initially  applied 
to  CPHF  theory  by  Handy  and  Schaefer  [41]  and  was  first  proposed  for 
many-body  derivative  equations  by  Adamowic,  et  al  [35].  It  was 
presented  systemmatically  by  Bartlett  [37]  and  Fitzgerald  et  al.  [36]. 
Salter,  et  al.  [42]  further  simplified  the  expression  to  a convenient 
computational  form. 

Recalling  that  1 = IPXPI  = iOXOj  + |$X#|,  the  derivative  energy 
can  be  written  as; 

AE^  = <0|(hJ  e^)^|0>  + <0lA  (hJJ  e^)^|0>  , (15) 

and  the  A equations  can  be  written  as 

_ fT» 

<0|A  (Hj^  - AE  <0|A|f>  + <0|(Hjj  e )^|$>  = 0 . 

Furthermore, 

<0|(Hj^  e*^)^  A|*>  = <0|(Hj^  e’’^)^|PXP|A|$> 
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= <0|(Hj^  e‘^)^|OXO|A|f>  + <Ol(Hj^  e’^)^i$Xf|A|f> 

= AE  <0|A|$>  + <0|(Hj^  $X$|  A|  ♦>  , 

so  the  the  A equations  may  be  written  in  terms  of  the  commutator 
[A, (Hj^e^)^] . As  in  the  case  of  the  coupled-cluster  equations,  only 
the  connected  product  is  non-zero.  Finally,  the  A equations  become; 

<0|(A  + <0|(Hj^e^)^|$X#|A|$>  + <0 1 (H^e^)^  | $>  = 0 . (16) 

Equations  (15)  and  (16)  are  the  final  expressions  for  the  the  CC 

y 

energy  derivative.  It  should  be  noted  that  AE^  is  expressed  in  terms 
of  unperturbed  t-amplitudes,  the  derivative  integrals  of  and  the  X- 
amplitudes.  Equation  (16)  may  be  used  to  determine  the  coupled  linear 
X-amplitudes.  The  equations  are  perfectly  general  and  are  not 
restricted  to  any  particular  model  of  coupled-cluster  theory  or  any 
choice  of  orbitals. 

The  general  equations  have  been  derived  and  presented  elsewhere 
for  the  full  CCSDT  model.  Here,  a canonical  Hartree-Fock  reference 
function  is  employed  for  computational  convenience.  The  remainder  of 
this  work  will  be  concerned  with  the  derivative  expressions  for  the 
coupled-cluster  doubles  (CCD)  model  and  MBPT  models  through  full 
fourth  order  (SDTQ-MBPT(4)) . As  stated  earlier,  low  iterations  of  the 
unperturbed  coupled-cluster  equations  provide  MBPT.  Additionally,  the 
low  order  iterations  of  the  coupled-cluster  derivative  equations 
generate  MBPT  derivatives.  Therefore,  in  this  work,  coupled-cluster 
theory  will  be  used  to  aid  in  the  identification  of  the  MBPT 
derivative  terms.  There  are  simplifications  that  occur  for  MBPT 


36 


derivatives  which  will  be  discussed  in  detail.  The  CCD  model  contains 
all  terms  required  to  evaluate  MBPT(2),  MBPT(3)  and  DQ-MBPT(4).  The 
SDTQ-MBPT(4)  model  is  complete  when  one  considers  the  case  of  the 
coupled-cluster  doubles  model  plus  linear  contributions  arising  from 
singles  and  triples  (CCD  + LST);  the  exponential  operator  is  truncated 
to  {exp  (T2)  + Tj  + T^}  for  this  approximation. 

The  operator  and  the  A operators  are  given  diagrammatically  in 
Figures  (7)  and  (8),  respectively.  The  fully-contracted  products 
which  contribute  to  the  CCD  and  SDTQ-MBPT(4)  derivative  energy 
expressions  are  represented  diagrammatically  in  Figure  (9).  The 
explicit  algebraic  equations  are  given  in  Table  (1).  Note  that  the 
algebraic  expressions  correspond  with  a given  diagram  and  are  so 
ordered. 

In  order  to  generate  all  terms  need  for  CCD  and  MBPT(4),  Equation 
(16)  represents  three  coupled,  linear  equations  for  the  A2,  and  A^ 
operators: 


<0|(A  (HNe'^)^)JS>  + <0|(Hj^e’^)JS>  = 0 

<0|(A  (Hjje'^)^,)^|D>  + <0|(Hjje^)^|SXS|A|D>  + <0 1 ( ) ^ | D>  = 0 

<0|(A  (Hjje'^)^)^|T>  + <0|(Hj^e'^)^|SXS|A|T>  + <0|  (Hj^e'^)^|DXD|  A|T>  = 0 . 

The  necessary  A equations  for  MBPT(4)  and  CCD  are  presented 
diagrammatically  in  Figures  (10-12).  The  corresponding  algebraic 
expressions  are  given  in  Tables  (2-4).  The  permutation  symbol 
(-l)P(pq)  means  the  identity  for  the  expression  plus  the  permutation 
of  the  labels  p and  q with  a corresponding  sign  flip  for  each 

P 

permutation  P.  A permutation  of  the  form  (-1)  P(p/qr)  amounts  to 
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three  terms  given  by  the  identity,  the  permutations  of  p with  q and  of 
p with  r.  The  expressions  are,  of  course,  a subset  of  those  given  by 
Salter  et  al.  [40]  for  CCSDT  and,  if  restricted  Hartree-Fock  (RHF) 
orbitals  are  chosen,  those  of  Scureria  and  Schaefer  for  CCSDT-1  [23]. 
As  in  the  coupled-cluster  energy  expression,  the  diagrams  contain  no 
closed  or  disconnected  parts. 

If  a comparison  of  the  X-amplitude  equations  with  the  unperturbed 
t-amplltudes  is  made,  it  becomes  apparent  that  several  X-amplitude 
equations  are  the  complex  conjugate  ("upside  down")  of  t-amplitude 
diagrams;  furthermore,  the  corresponding  algebraic  expressions  are 
identical  for  real  orbitals.  In  fact,  the  amplitudes  are  identical 
for  any  linear  model  of  many-body  theory  (MBPT(2),  MBPT(3),  SDT- 
MBPT(4),  LCCD,  etc.).  Since  the  equations  are  linear,  the  X- 
amplitudes  and  t-amplitudes  differ  due  to  the  the  quadratic  terms  of 
T2  (diagrams  5a-d  in  Figure  (5)).  The  X-amplitudes  contain  linear, 
yet  similar  terms  (diagrams  11-17)  in  Figure  (11).  If  the  quadratic 
contributions  to  T2  are  given  by  Q®^(I,txt),  the  A2  contributions  may 
be  expressed  Q®J(t,IxX)  + Q®^(t,XxI).  This  relationship  has  great 
consequence  to  MBPT  because  only  the  quadruples  of  MBPT(4)  require 
special  treatment. 

For  the  fourth-order  singles,  the  linear  contributions  of  to 
K^,  diagrams  (18,19)  in  Figure  (11),  are  analogous  to  the  linear 
contributions  of  T^^  to  T2,  represented  by  diagrams  (3a, b)  in  Figure 
(5).  For  fourth-order  triples,  the  linear  contribution  of  A^  to  A2, 
diagrams  (20,21)  of  Figure  (11),  are  analogous  to  diagrams  (4a-b)  of 
Figure  (5).  The  linear  contributions  of  A^  and  to  Aj^  are  given  by 
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diagrams  (1,2)  and  (3,4)  in  Figure  (10).  Likewise,  the  corresponding 
and  T2  contributions  to  are  given  by  diagrams  (3a, b)  and  (2b, c) 
of  Figure  (4).  Finally,  there  is  the  linear  contribution  of  and 
to  which  is  represented  by  diagrams  (24,25)  and  (22,23)  in  Figure 
(12);  Analogously,  the  T2  and  linear  contributions  to  are  given 
by  diagrams  (la,b)  and  (3a, b)  of  Figure  (6). 


Elimination  of  Singularities 

To  summarize  to  this  point,  the  general  derivative  expressions 

have  been  given  for  CCD  and  MBPT(4).  The  need  for  3N  derivative  t- 

amplitudes  has  been  eliminated.  Also,  the  relationship  of  the 

perturbation  independent  X-amplitudes  to  the  unperturbed  t-amplitudes 

has  been  discussed.  At  this  point,  it  is  necessary  to  discuss  the 

remaining  perturbation-dependent  terms  arising  from  the  derivative 

y y 

integrals  in  the  f^  and  operators. 

It  is  appropriate  to  examine  a particular  orbital  basis  function 
♦p(X)  contained  in  the  set  {<|'(X)}.  The  basis  function  is  normally  a 
linear  combination  of  primitive  atomic  orbital  functions  {<|'(X)}  and 
expansion  coefficients  C (MO  coefficients  obtained  from  Hartree-Fock 
theory) ; 

♦p(X)  = |p>  = l'KX)C(X)  . 

The  derivative  of  the  above  expression 


♦p«<) 

8X 


3'KX) 

ac(x) 

s 

> 

C(0)  + 

|«|'(0)>  > 

x*o 

ax 

x=o 

ax 
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y { p^q}  = ^ 

<pq||rs>’‘{p+qV}  = 

. J„|  . ] 

. {..y  . j...|  . J..J 

y 

Figure  7.  The  Diagrammatic  Representation  of 


K'b{'Vb}= 

x'i^jitarbktc}=  /gg\ 


% * A'A 

* f"A 


Figure  8.  The  Diagrammatic  Representation  of  the  De-excitation 
Operators  ^2,  and  A^. 
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Figure  9.  The  Diagrammatic  Representation  of  the  Coupled-Cluster 
Energy  Derivative. 
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A-  A-  IP  IP 

I 2 3 4 

Figure  10.  The  Diagrammatic  Representation  of  the  Equation. 
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Figure  11.  The  Diagrammatic  Representation  of  the  Equation. 
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Figure  12.  The  Diagrammatic  Representation  of  the  Equation. 
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Table  1. 
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The  Derivative  Energy  Expression 


1 

5 

ab 

<ij| |ab>^ 

x'  f 

a 

a 

j 

x'  f 

a 

ab'^' 

b 

i 

1 

2 

X^  <ai| |bc>^  t 

be 

ji 

1 

2 

xj  <jk|  |ia>^  t 

ba 

jk 

1 

5 

ab 

<ab| |ij>^ 

1 

2 

<ab 1 1 ci>^ 

1 

2 

ab 

<ka| |ij>^ 

1 

f X 

ab 

2 

ab 

'kj 

1 

2 

x^j 

ac 

ab 

'5S 

1 

8 

ab 

<ab| 

cd>^ 

1 

8 

ab 

<kl| 

ij>^ 

'u 

ac 

<ja| 

bi>^ 

'w- 

kj 

1 

? 

<ai| 

bc>^ 

dbc 

^kji 

1 

5 

<jk| 

ia>X 

eba 

^Ijk 

1 

2 

cd 

<ij| 

ab>X 

ki 

db 

Ij 

1 

5 

cd 

<ij| 

ab>X 

t?® 

ij 

db 

kl 

1 

16 

cd 

<ij| 

ab>X 

^kl 

cd 

ij 

1 

5 

cd 

<ij| 

ab>X 

ki 

1 

5 

x‘y  <ab 
abd 

|ci> 

X .dc 
^kj 

1 

5 

X^^L  <ka 
acb 

|ij> 

X .cb 
hk 

for  CCD  and  MBPT(4). 
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Table  1 — continued 


1 

.ilk 

f X .abc 
*ji  hkj 

■ 12 

abc 

1 

^ 12 

xljk 

adc 

f X .deb 
ab  ^ijk 
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Table 
0 = 


Note: 


2.  The  Equation  the  CCD  and  MBPT(4), 


ba 


f.. 
a Ij 


xJJ  <bc| |aj> 
Xja  <ib||jk> 


1 

2 

1 

2 

Summation  over  all  but  the  target  labels 
i and  a is  implied. 
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Table  3.  The  Equation  for  CCD  and  MBPT(4). 
0 = <ij I I ab> 


+ 

P(lj) 

\ <cj 1 1 ab> 

- 

P(ab) 

<i j 1 1 

tb> 

+ 

P(ij) 

f-i, 
ab  jk 

- 

P(ab) 

f . 
ca  cb 

+ 

1 

2 

j <cd 1 
cd  ' 

1 ab> 

+ 

1 

2 

^ab  1 

lkl> 

+ 

P(ab|ij) 

<jc| 
ca  ' 

|bk> 

- 

^ p(ij) 

\ki  . 1 
\d  1 

1 UK  ..cd 
|ab>  t^^ 

+ 

1 

5 

*cd 

|ab>  t^5 

- 

2 P(ab) 

'‘ca  1 

idb> 

+ 

P(ab|ij) 

X*^^  <lj| 
ca  ' 

idb> 

- 

2 P(ab) 

X^^  <kl| 
ca  ' 

|db>  t^J 

+ 

1 

1 

,kl  . 1 
\b  1 

1 jv  »cd 

ici> 

- 

J p(ij) 

\ki  >1  . 1 
\b  1 

1 JK  ..cd 

|cd> 

+ 

2 P(ab) 

<cd 

cda 

1 |bk> 

- 

2 P(ij) 

,kli 

\ab 

1 |kl> 

Note: 

Summation  over  all  but 
and  b is  implied. 

the  target 

P(pq)  is  an  operator  that  generates  two  terms: 

P(pq)  = identity  - interchange  p«»q. 

P(ab|ij)  is  an  operator  that  generates  four  terms: 
P(ab|ij)  = identity  - interchange  a«b  - interchange  i«j 
+ interchange  both  a«b  and  i<>j. 
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Table  4. 
0 * - 
+ 


+ 

Notes  Summation  over  all  but  the  target  labels  i,  j,  k,  a,  b and 
c is  implied. 

P(p/qr)  is  an  operator  that  generates  three  terms: 

P(p/qr)  =»  identity  - interchange  p«>q  - interchange  p^r. 

P(a/bc|l/jk)  is  an  operator  that  generates  nine  terms: 
P(a/bc|i/jk)  = identity  - interchange  a®b  - interchange  a»c 
- interchange  i«j  - interchange  i<»k 
+ interchange  both  a»b  and  i<»j  + interchange  both  a<»b  and  i«»k 
+ interchange  both  a<»c  and  i®j  + interchange  both  a«>c  and  i«k. 


The  Equation  for  CCD  and  MBPT(4). 
P(a/bc|k/ij)  <dk| |bc> 
P(c/ab|i/jk)  Xll  <jk| |lc> 

'kl 

P(c/ab) 
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contains  two  parts.  The  first  is  the  derivative  of  the  primitive 
functions  (AO  part)  multiplied  by  the  known  unperturbed  MO 
coefficients;  the  second  is  the  derivative  of  the  SCF  molecular 
orbital  coefficients  (MO  part)  multiplied  by  the  primitive  functions. 
The  MO  part  is  responsible  for  orbital  rotation  and  may  be  obtained 
from  solutions  of  Coupled-Perturbed  Hartree-Fock  (CPHF)  theory  [2,43]. 

According  to  Gerratt  and  Mills  [2],  the  first-order  change  of  a 
molecular  orbital  subject  to  a perturbation  is  given  by: 


♦p(x)  . +p(0)  * £ ujp  •t'p«» 

or,  specifically, 


d|a> 


ax 

X=0 


X 

|a  > 


XX  X 

EU  |f>+  SU  |k>+  EC  |u> 
fa  ka  ya 

tta  k y 


(17) 


3|i> 


ax 

x-o 


X 

11  > 


XX  X 

EU|f>+EU|k>+EC|y> 
fi  ki  yi 

f k/i  y 


(18) 


where  the  labels  y,v,X, o, ...  are  used  to  represent  the  primitive  AO 
basis  functions.  The  MO  coefficients  for  the  unperturbed  orbital  p 
are  given  by  C^p.  The  coefficients  are  the  CPHF  coefficients  and 
are  subject  to  the  orthogonality  condition: 


X X * X 

U + (U  ) + S = 0 . (19) 

pq  qp  pq 

As  defined  by  the  derivative  of  the  SCF  orthonormality  condition  <p|q> 
= 6p^.  Simply  stated,  the  derivative  of  an  SCF  orbital  at  X=0,  is 
given  by  a linear  combination  of  the  unperturbed  SCF  orbitals  and  a 
term  arising  from  the  derivative  of  the  primitive  AO  basis  functions. 
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The  CPHF  equations  consists  of  three  distinct  blocks:  (i)  the 
virtual-occupied  block,  (ii)  the  occupied-occupied  block,  and  (iii) 
the  virtual-virtual  block.  The  virtual-occupied  block  of  the  CPHF 
coefficients  are  determined  by  solving  the  linear  equation: 


X 

Q 

X -1  bj 

U = E (A  ) 

ai  ai,bj  e - e 

bj  j b 


(20) 


where 


<ab| I ij>  + <aj I I ib> 

A = 1 + 

ai,bj  e - e 

a i 


and  Ep  is  the  unperturbed  orbital  eigenvalue  of  orbital  p.  Given  the 
virtual-occupied  block,  the  occupied-occupied  and  virtual-virtual  blocks 
are  then  determined: 


X -1  X X X * 

U = [ Q + E { U <km||ig>  + (U  ) <kg||im>}  ] (21) 

ki  (e  -e  ) ki  gm  gm 

k i g,m 


X -1  X X X * 

U = [ Q + Z { U <fm|  |ag>  + (U  ) <fg|  |am>}  ] 

fa  (c  -E  ) fa  gm  gm 

f a g,m 

the  diagonal  elements  of  the  Fock  matrix  are  given  by: 

XX  X X * 

E = Q + Z { U <am| |ag>  + (U  ) <ag| |am>} 
a aa  gm  gm 

g,m 


(22) 


XX  X X * 

E = Q + Z { U <im||ig>  + (U  ) <ig||im>) 
i ii  gm  gm 

g.m 


Note  that 
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V 

The  elements  of  the  general  one-electron  perturbation  matrix  Q are 
given  by: 


XXX  X 

Q - h - S C - Z S <pl|  |qk> 

pq  pq  pq  pq  ^i 

k,i 


* 3<uX|  I V(T> 

+ E (C  ) C P 

UP  N)q  Xcr  3X 
U , V X=0 

X,  a 


(23) 


where 


c 

pq 


' 1 when  p and  q are  both 

- (e  + e ) 

2 p q virtual  or  both  occupied 


e otherwise 

q 


In  the  above  equation,  represents  a matrix  element  of  the 

y 

Hartree-Fock  density  matrix  in  the  AO  basis.  S is  the  derivative  of 
the  overlap  matrix: 


X * 3<u  I v> 

S = E (C  ) C 

pq  up  vq  3X 

U,v  x=o 


and  h^ 


is  the  derivative  of  the  core  Hamiltonian  matrix: 
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X 


9<y i h 1 v> 


h 


pq 


z (C  ) c 

yp  vq  3X 

U,v 


X=0 


Equation  (22)  is  a generalization  of  equation  (51)  in  Reference  [5]. 

Early  derivative  studies  employed  the  above  equations  to  formulate 
a final  expression  for  the  derivative;  however,  the  use  of  the 
occupied-occupied  and  virtual-virtual  CPHF  coefficients  lead  to 
singularities  for  systems  that  contain  degenerate  or  near-degenerate 
orbitals  due  to  the  form  of  the  denominator  (£^  - 

As  pointed  out  by  Handy,  et  al.  [44],  "In  large  basis  set  calculations 
or  in  molecules  with  degenerate  symmetry  groups  this  problem  will 
invariably  occur."  However,  this  problem  can  be  avoided  by 
factorization  of  the  derivative  expression  or  by  choice  of  non- 
canonical  Hartree-Fock  orbitals.  In  this  work,  the  latter  choice  is 
incorporated  and  the  factorization  approach  is  given  elsewhere  [42]. 

The  idea  of  non-canonical  orbitals  in  MBPT(2)  gradients  was  first 
employed  by  Handy;  although,  the  fundamental  approach  was  found  in 
early  derivative  studies  of  Moccia  [45]  and  Pulay  [2].  In  this 
approach  the  occupied-occupied  and  virtual-virtual  CPHF  coefficients 

are  given  by: 


U 


X 


1 X 
- S lf> 


(24) 


fa 


2 fa 


U 


X 


ki 


1 X 

- S I k>  , 

2 ki 


(25) 
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This  approach  is  valid  since  CC  and  MBPT  are  invariant  to  orthogonal 

transformations  of  the  occupied  or  virtual  orbitals.  Now  the  use  of 

the  pathological  denominators  is  avoided  in  the  derivative  expression. 

For  the  choice  of  non-canonical  orbitals,  however,  the  off-diagonal 

elements  f . ^ and  f . are  not  zero.  These  matrix  elements  are  given 
ab  1 j 

by; 

f «=  + Z (<am||bg>  + <ag||bm>)  (26) 

ab  ab  gm 

f.,^  * 0^,  + Z uj|^(<im|  |jg>  + <ig|  I jm>)  . (27) 

^ J g,m  ® 

The  expressions  in  Equation  (19)  and  (24-27)  may  now  be  applied  to 
the  derivative  expression  in  Table  (1).  For  convenience,  the 
derivative  energy  expressions  may  be  grouped  together  based  on  the 
common  derivative  integrals  and  derivative  Fock  matrix  elements: 


* [j^lj  ^u^ab  ^ab^ 


ab 


+ Z F(ij,ab)  <ij||ab>^  + Z F(ab,cd)  <ab||cd> 

i,j,a,b  a,b,c,d 

+ Z F(kl,ij)  <kl| |ij>^  + Z F(ja,bi)  <ja| |bi>^ 

i,j,k,l  i,j,a,b 

+ Z F(ai,bc)  <ai| |bc>^  + Z F( jk, ia)  <jk| | ia>^  , 

i,a,b,c  i>j>k.)a 


D and  T intermediates  for  the  CCD  and  SDTQ-MBPT(4)  models  are  defined 
in  Tables  (5)  and  (6),  respectively.  It  should  be  noted  that  F(ai,bc) 
and  F(jk,ia)  are  identically  zero  for  approximations  limited  to  T2 
(i.e.  double  excitations). 
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The  derivative  energy  expression  is  expanded  in  terms  of  the  CPHF 
identities. 

AE^  = (Qjj  + (<ie||jm>  + <im||je>)) 

+ D . (q\  + (<ae||bm>  + <am|lbe>)) 

ab  ' ab  em  ' ' ' 

+ (r(ij,ab)  - r(ji,ab))  (U^.<ej|lab>  - 1/2  sJJ. <mj  | | ab>) 

+ (r(ij,ab)  - r(ij,ba))  S^^<ij||eb» 

+ (r(ab,cd)  - r(ba,cd)  + r(cd,ab)  - r(cd,ba)) 

X ([-U^  -S^  ]<mb||cd>  - 1/2  <eb||cd>) 

am  am  ea 

+ (r(kl,ij)  - r(lk,ij)  + r(ij,kl)  - r(ij,lk)) 

X (uX^<el||ij>  - 1/2  S^^<mll|ij» 

+ (r(ja,bi)  + r(ib,aj))  (U^j<ea||bi>  - 1/2  S^Xma||bi>) 

+ (r(ja,bi)  + r(ib,aj))  (l-U^m-San,]<jni||bi>  - 1/2  S^^<je||bi» 
+ r(ai,bc)  (l-U3^-Sa„,l<'"i|  |bc>  - 1/2  S^^<ei||bc>) 

+ r(ai,bc)  (U^.<ae||bc>  - 1/2  S^.<aml |bc>) 

+ (r(ai,bc)  - r(ai,cb))  sjj^<ai  | |ec>) 

+ (r(jk.ia)  - r(kj,ia))  (U\<ek||ia>  -1/2  s\<mk|lia>) 

+ r(jk,ia)  (u\<jk||ea>  - 1/2  S^.<jk|  lma>) 

+ r(jk,ia)  ([-uj  -S^^l<jk||im>  - 1/2  <jk||ie» 


+ AO  derivative  part 
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The  contributions  arising  from  the  derivative  of  the  primitive  AO 
basis  functions  are  not  explicitly  shown  here  because  the  remainder  o± 
the  discussion  concerns  the  remaining  terms.  The  final  form  of  the  AO 
derivative  part  is  given  in  Table  (7).  Following  Fitzgerald  at  al. 
[35],  r(ywpff)  represents  the  F intermediates  transformed  to  the 
primitive  basis  set.  This  avoids  the  transformation  of  3N  set  of 
derivative  integrals  to  the  molecular  orbital  basis.  A dot-product 
between  the  r(uvp<y)  intermediate  and  the  3N  derivative  integrals  may 
be  performed  as  the  derivative  integrals  are  produced  for  each  degree 
of  freedom. 

The  remaining  terms  or  the  MO  derivative  part  then  takes  the  the 
following  form  where  only  the  virtual-occupied  block  of  the  CPHF 
coefficients  remain. 


MO  derivative 

- E D 

+ E D + 

2 E X 

M 'J 

. ab  ab 

ai  ai 
a,i 

part 

a,b 

* ^ij 

* \ lab 

2 lai  • <28) 

a,i 

i,j  ^ 

a,  b 

Recalling  the  form  of  the  CPHF  coefficients,  the  Z-vector  method 
may  applied  to  the  third  term  of  Equation  (28).  This  term  may  now  be 
expressed  in  equivalent  form: 


X 

Q 

X -1  bj 

2EX  U = 2EX  E(A)  

ai  al  ai  ai,bj  t -t 

a,i  ai  bj  j b 


X 

2 E D Q 
bj  bj 
b,  j 


(29) 


where  the  virtual-occupied  block  of  the  matrix  is  the  solution  to 
the  perturbation-independent  linear  equation: 
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ZD(e-e)A  =X  . (30) 

bj  j b bj,ai  ai 

bj 


Therefore  the  MO  derivative  part  becomes: 


XX  X 

MO  derivative  « I D Q + Z D Q + 2ZD  Q (31) 

ij  ij  ab  ab  ai  ai 

part  i,j  a,b  a,i 


XXX 
+ EIS  + ZIS+2ZIS 

ij  ij  ab  ab  ai  ai 

i,j  a,b  a,i 

Equation  (31)  is  in  final  computational  form. 

The  finite-order  MBPT  derivatives  may  be  made  explicit  from  the 
above  expressions.  In  order  to  compute  a finite-order  derivative  of 
order  n in  correlation  energy,  the  n^^-order  MBPT  density  matrix  D and 
the  n''^-order  I intermediates  are  required.  The  algebraic  formulas 
for  MBPT(2),  MBPT(3),  MBPT(4),  D-MBPT(”)  (CCD  restricted  to  linear 
terms  or  LCCD)  and  CCD  are  given  in  Tables  (9-13).  After  appropriate 
construction  of  the  X,  the  virtual-occupied  block  of  the  density 
matrix  is  obtained  by  solution  of  Equation  (30).  The  expressions  are 
given  in  terms  of  the  low-order  iterates  of  T and  the  function 
Q(A,BxC).  The  equations  for  the  low-order  iterations  of  A and  T are 
given  in  Table  (14). 

In  the  CCD  expression  there  are  four  terms  of  the  form  Q(A,BxC): 
(1)  the  quadratic  contribution  to  the  unperturbed  energy,  (2)  Q(t,Xxt) 
and  (3)  Q(t,IxX)  which  together  contribute  to  A,  and  (4)  Q(X, txt) 
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which  is  defined  in  the  r(ij,ab)  expression.  Due  to  the  fact  that  A = 
T for  low-orders  of  perturbation  theory,  these  four  terms  are  closely 
related.  It  is  shown  in  the  Appendix  that  only  two  of  these  terms  are 
required  for  the  computation  of  the  quadruples  contribution  of  MBPT 
because  the  following  relationship  may  be  derived: 

Q(I,txt)  + Q(t,txt)  = Q(t,Ixt)  + Q(t,txl)  . 


Identification  of  the  Many-Body  Response  Density 


It  is  well  known  that  the  construction  of  a density  matrix  [45] 
from  an  approximate  electronic  wavefunction  is  a convenient  way  to 
treat  many  one-electron  properties.  The  one-electron  reduced  density 
matrix  is  defined  as 


'f(l...n)'f'*(l'...n)dT2...dTj^ 


* 

= E D ♦ (1)  (1')  • 

p,q  pq  p q 


(32) 


The  electron  density  is  the  "diagonal"  of  the  reduced  density  matrix 

•k 

y(1,1')|  = p(l)  = i:  D ♦ (D-t-  (1)  . 

1-1'  p,q  pq  p q 

The  expectation  value  of  some  one-electron  operator,  Q(l)  is  given  by 


<Q> 


[Q(l)p(l|l')l  dT  . 

1=1'  1 
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Table  5.  The  Occupied-Occupied  and  Virtual-Virtual  Blocks  of  the 
CCD  and  MBPT(4)  Density  Matrix. 


Dij  . - 

tf 

a i 

1 

2 

Jk  ab 
ab  ’^ki 

1 

.jlk  abc 

■ 12 

\bc  ^Iki 

D w = 
ab 

t*’ 
a 4 

1 

■ 2 

ac  ij 

1 

12 

.ijk  deb 
adc  Mjk 

Note:  Summation  over  all  but  target  labels  is  implied. 
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Table  6.  The  F Intermediates. 


F(ij ,ab) 

s 

X 

5 

^ i 

ab 

1 

ab 

+ 

5 

'ij 

1 

ca 

.db 

+ 

2 

\d 

^ki 

'ij 

1 

.ca 

,db 

? 

^d 

t . . 
ij 

^kl 

1 

ab 

cd 

+ 

15 

^d 

^kl 

'ij 

1 

ab 

cd 

? 

^d 

^ki 

'ij 

F(ab,cd) 

1 

8 

ab 

F(kl,ij) 

_ 

1 

8 

ab 

F(ja,bi) 

- 

ac  kj 

F(ai , be) 

= 

1 

2 

a Ji 

- 

1 

5 

.kJ  dbc 
ad  ’^kji 

- 

1 

2 

be  j 

+ 

1 

? 

xikj  da 
%cd 

F(jk,ia) 

= - 

1 

2 

+ 

1 

5 

xil  eba 
cb  hjk 

+ 

1 

2 

ab  1 

- 

1 

5 

Jkl  cb 
acb  hi 

Note:  Summation  over  all  but  the  target  labels  is  implied. 
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Table  7.  The  AO  Derivative  Part  of  the  Coupled-Cluster  and  MBPT 
Energy  Derivatives. 


AO  derivative  part  » 2)  <w(X)v(X)  | p(X)  r(uvpo)  , 

y,v, 

p,  a 


where 


rcuvp,)  . I r(ij,ab)  (c^  c*j  c^. 

i » J 

a,b 

. £ r(ab,cd)  4 

a I D 

c,d 


. £ r(ij.ki)  c,j  - c,,^)  c^j 

i f J 

k,l 

. E r(ja,bi)  (Cp^  C^,  - Cp.  C^J,)  c^.  c;^ 
i,a 
j>b 

. E r(ai,bc)  (C  C^^  - C C^^)  c;^  C^. 

i,a 


b,c 


★ ★ 


. £ r(jk,ia)  (C^J  C^j)  C^. 

i , a 
j»k 


Note:  <yv|  p(T>^  denotes  the  derivative  of  an  ordinary  AO  two-electron 

integral,  evaluated  at  X=0. 
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Table  8.  The  I and  X Intermediates  Used  in  the  Evaluation  of  CCD 
and  MBPT  Energy  Derivatives. 


- -1/2  (r(jk,ab)  - r(kj,ab))  <ik||ab> 

-1/2  (r(jl,mk)  - r(lj,mk)  + r(mk,jl)  - r(mk,lj))  <il||mk> 
-1/2  (r(ja,bk)  + r(kb,aj))  <ia| |bk> 

-1/2  r(aj,bc)  <ai| |bc> 

-1/2  (T(jk,la)  - r(kj,la))  <ik||la> 

-1/2  r(lk,ja)  <lk| |ia> 

. -1/2  (r(ij,bc)  - r(ij,cb))  <ij| |ac> 

-1/2  (r(be,cd)  - r(eb,cd)  + r(cd,be)  - r(cd,eb))  <ae||cd> 
-1/2  (r(jb,ci)  + r(ic,bj))  <ja||ci> 

-1/2  r(bi,cd)  <ai| |cd> 

-1/2  (F(di,bc)  - r(di,cb))  <di||ac> 

-1/2  r(jk,ib)  <jkl |ia> 

I . = -1/2  (r(kj,ab)  - r(kj,ba))  <kj|lib> 

ai 

-1/2  (r(ab,cd)  - r(ba,cd)  + r(cd,ab)  - r(cd,ba))  <ib|lcd> 
-1/2  (T(ja,bk)  + T(kb,aj))  <ji||bk> 

-1/2  r(aj,bc)  <ij||bc> 

-1/2  (r(dj,ac)  - r(dj,ca))  <dj||ic> 

-1/2  r(jk,la)  <jk|lli> 
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Table  8 — continued 


+1/2  (<ka||ji>  + <ki||ja>) 

+1/2  (<ca||bi>  + <ci||ba>) 

+1/2  (r(ij,cb)  - r(ji,cb))  <ajl|cb> 

+1/2  (r(il,kj)  - r(li,kj)  + r(kj,ll)  - r(kj,li))  <al||kj> 
+1/2  (r(ic,bj)  + r(jb,ci))  <ac||bj> 

+1/2  r(di,bc)  <da| |bc> 

+1/2  (r(ik,jb)  - r(ki,jb))  <ak| | jb> 

+1/2  r(jk,ib)  <Jk| |ab> 


Note;  Summation  over  all  but  the  target  labels  is  implied. 
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Table  9.  Summary  for  the  MBPT(2)  Model. 


D.  . = 

IJ 

1 

2 

D u •=  - 

ab 

1 

2 

tjjd)  tjJ(l) 

r(ij,ab) 

s 

1 

2 

all 

other  F's  = 0 

Note:  Summation  over  all  but 

target  labels  is  implied. 

Table  10. 

Summary 

for  the  MBPT(3) 

Model. 

“ij  ■ 

1 

2 

t^^l)  tjJ(2) 

1 

^ 2 

°ab  = - 

1 

2 

t®J(l)  tjJ(2) 

1 

" 2 

At^J(2)  tjjd) 

r(ij,ab) 

1 

2 

r(ab,cd) 

= 

1 

8 

t^j(l)  t^^ 

(1) 

r(ki,ij) 

= 

1 

8 

(1) 

r(ja,bi) 

= 

(1) 

; r(ai,bc)  = r(jk,; 

Note:  Summation  over  all  but  target  labels  is  implied. 
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Table  11.  Summary  for  the  MBPT(4)  Model. 

= - tj(2)  t®(2) 

* 2 'k^3)  * 2 1 t'J(l) 

+ ^ Qj^(t,  t(l)xt(D)  t^^(l) 

- 12 

Dab  = '?<') 

- 2 '?j(»  'u»)  - 2 - 2 

- I Q®J(t,  t(l)xt(D)  t^J(l) 

12 


r(ij,ab)  = 

1 

2 

,ab,~.  1 

tij(3)  + 2 

Q^^(td),  t(l)xtd)) 

r(ab,cd)  « 

1 

8 

tjjd)  tj^(2) 

1 

8 

M^J(2)  tj^(l) 

r(kl,ij)  = 

1 

8 

tjjd)  t®J(2) 

1 

^ 8 

At^j(2)  t®J(l) 

r(ja,bi)  = 

tllm  t'5(2) 

+ 

At-(2)  t^^b(l) 

r(ai,bc)  = 

t®(2)  t^Jd) 

+ 

2 tj^d) 

r(jk,ia)  = 

- 

tj(2) 

- 

2 tjjd) 

Note:  Summation  over  all  but  target  labels  is  implied. 


63 


Table  12.  Summary  for  the  LCCD  (or  D-MBPT(®))  Model. 


Dij  = 

1 

2 

ab. 

-) 

\b  = - 

1 

2 

t®'^( 

”)  t'Jc 

>) 

r(lj,ab) 

= 

1 

2 

r(ab,cd) 

s 

1 

8 

tj5(«) 

r(ki,ij) 

= 

1 

8 

tfj(-) 

r(ja,bi) 

3: 

1k<"> 

t^J(»)  ; 

r(ai,bc)  = r(jk,ia) 

Note:  Summation  over  all 

but  target 

labels  is  implied. 

The  symbol  represents  the  converged  LCCD  amplitudes. 
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Table  13. 

Summary 

for 

the 

CCD 

°ab  = 

1 Jk 

2 \b 

‘^ab=  - 

2 ac 

'u 

r(ij ,ab) 

1 

? 

ab 

1 

* 4 

ab 

1 

2 

cd 

'u 

db 

1 

■ 5 

cd 

db 

^kl 

1 

^ 15 

ab 

’^kl 

'u 

1 

■ 5 

'u 

r(ab,cd) 

1 

8 

ab 

'u 

r(kl,ij) 

1 

8 

ab 

’^kl 

r(ja,bi) 

B 

ac 

cb 

^kj 

Note;  Summation  over  all  but  the  target  labels  is  implied. 


65 


Table  14.  The  Low-order  Iterates  of  A and  T. 


First  order: 


X^j(l)  = <ij||ab>/D^J 

nab 

= Ei^  e.- 

Second  order: 

tjj(2)  = 

X‘J(2)  . t»5<l)  . At=J(2) 

At®^(2) 

= [\  tj^(l)  <ab||cd>  + \ t®J(l)  <ij||kl> 

+ (-l)^P(abjij)  t^^(l)  <jcl|bk>  ] / 

Third  order: 

■ 

t®J(l)  + At®^(2)  + At®^(3;SDT)  + At^^(3;Q) 

^ab<3)  - 

t®J(l)  + At®J(2)  + At®^(3;SDT)  + AX^j(3;Q) 

At®J(3;Q)  - Q®J(I,  t(l)xt(D)  / 

iX*j(3,Q)  - 

{Q®^(t(D,  Ixt(D)  + Q?^(t(l),  t(l)xl)}  / uf. 

= 

At®^(3;Q)  + Qij(t(l),  t(l)xt(D) 

t(l)  = t®j(l)  1=  <ij||ab> 

Q®J(A,  BxC)  = 

1 cd  r -cd  „ab  „ f -ac  pbd  gbd  ac  'l 

5 \l  L ®ij  ^kl  - M ®ij  ^kl  ^ ®ij  ^kl  j 

- 2 1 

f „ab  _cd  -ab  „cd  'j  , f „ac  „bd  bd  _ac  'I  ] 

1 ®ik  'jl  * ®J1  J * 1 \i  '^Ij  * ®ki  '^Ij  > \ 
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Table  14 — continued 

Atjj(3;SDT)  . [ ^ Atjj(2)  <ab| |cd>  + ^ At®J(2)  <ij | | kl> 

+ (-l)^P(ablij)  At^i(2)  <jc||bk> 

+ (-l)^P(ij)  tj(2)  <cj||ab>  - (-l)^P(ab)  t®(2)  <ij||kb> 


+ 1 (-l)^P(ab)  t^i|(2)  <cd||bk> 

- 

\ (-l)^P(ab)  t^®J(2)  <jc||kl>  ] / 

tj(2) 

= X^(2)  - [ - ^ t^Jd)  <bc|  |aj>  + 

1 tj^(l)  <ib||jk>  ] / D® 

D?  = e.-  e 

1 1 a 

= X^J^(2)  - [ - (-l)^P(a/bclk/ij) 

tj^(l)  <dk||bc> 

+ (-l)^P(c/ab|i/jk)  tl\n)  <jk||lc>  ] / 

_abc 

D...  = e.+  e.+  - e - 

ijk  1 j K a b c 


Note:  Summation  over  all  but  target  labels  is  implied. 
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In  the  occupation  number  representation,  the  matrix  elements 

D = <^18  a \Y>  (33) 

pq  P q 

constitute  the  discrete  representation  of  the  one-electron  reduced 
density  matrix  while  equation  (32)  is  the  continuous  representation. 

The  expectation  value  approach  to  one-electron  properties  leads  to 
complications  in  many-body  methods.  For  example,  in  coupled-cluster 
theory,  the  exponential  ansatz  leads  to  a non-terminating  series  in 
Equation  (33).  In  MBPT  models,  one  can  define  a density  by  selecting 
a given  MBPT  vavefunction  or  by  truncating  the  perturbation  expansion 
of  the  vavefunction  to  a given  order  in  the  overall  correlation  (in 
this  case  there  is  no  corresponding  vavefunction  associated  with  the 
MBPT  model).  Using  an  SCF  reference  function,  would  consist  of 

only  double  excitations,  and  the  density  matrix  constructed  from  it 
would  not  give  meaningful  results  for  the  second-order  density,  since 
it  would  neglect  all  single  excitation  effects  which  are  crucial  to 
the  correct  description  of  one-electron  properties.  Therefore,  a 
better  approximation  to  the  second-order  density  would  be 

, (1)  + (1)  (0)  + (2) 

D"^=<Y  laalY  > + <T  laalf  > (34) 

MBPT  p q MBPT  MBPT  p q MBPT 

(2)  + (0) 

< Y I a a I Y > , 

MBPT  p q MBPT 

(2) 

which  includes  selected  single-excitation  contributions  from 
[46]. 

Since  the  Hellman-Feynmann  theorem  is  not  satisfied  for  many-body 
methods  (i.e.  the  energy  derivative  is  not  equal  to  the  expectation 
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value),  one  can  define  a given  property  as  either  an  expectation  value 
or  an  energy  derivative.  The  energy  derivative  approach  allows  for  an 

I 

unambiguous  definition  of  a density  matrix  and  provides  a closed-form 
expression  for  the  density  matrix.  In  this  work,  the  density  matrix 
defined  in  this  manner  is  referred  to  as  a "response"  density. 

In  addition  to  actual  correlation  corrections,  the  response 
density  includes  all  effects  due  to  orbital  relaxation  that  arise  from 
orbitals  changing  to  accommodate  the  perturbation;  these  effects  are 
included  via  the  coupled-perturbed  Hartree-Fock  coefficients. 
Relaxation  and  actual  correlation  effects  are  separable  in  a 
diagrammatic  analysis  and  have  been  characterized  [41,47,48].  The 
first  two  papers  focused  on  second-order  properties;  the  effect  of 
orbital  relaxation  in  the  CCSD  model  for  first-order  properties  was 
discussed  by  Salter  et  al.  [48].  Furthermore,  it  has  been  shown  that 
such  relaxation  effects  are  partly  responsible  for  the  satisfaction  of 
the  Hellman-Feynmann  theorem  for  Cl  wavefunctions  [49],  which  in  that 
case  manifest  themselves  in  terms  of  introducing  Cl  triple 
excitations. 

An  effective  density,  which  is  independent  of  any  perturbation, 
was  introduced  in  CC  gradient  theory  [37,38].  The  concept  has  been 
generalized  to  define  a response  density  by  explicitly  separating  the 
AO  contributions  from  the  MO  contributions. 

Traditional  many-body  methods  for  evaluating  one-electron 
properties  have  employed  finite-field  techniques.  Although  reliable, 
these  techniques  rely  on  the  numerical  construction  of  the  energy 
derivative  with  respect  to  an  external  field;  this  normally  requires 
at  least  two  field  values  for  each  property  of  interest.  A much  more 
efficient  route  should  require  only  a single  correlated  calculation. 
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A direct  method  utilizing  equation  (33)  has  been  employed  for  some 
MBPT  examples  [50].  However,  as  discussed  above,  there  are 
ambiguities  in  this  approach  since  the  neglect  of  orbital  relaxation 
can  be  important.  The  response  density  method  outlined  below  requires 
a single  correlated  calculation  and  includes  all  orbital  relaxation. 

Response  properties  are  associated  with  the  derivative  of  the 
energy  with  respect  to  a perturbation  parameter  X*  This  parameter 
controls  the  magnitude  of  a one-electron  potential  energy  perturbation 
XQ.  For  example,  the  z-component  of  the  dipole  moment  is  obtained  by 
introducing  a uniform  external  electric  field  in  the  z-direction,  and 
Equation  (23)  reduces  to; 


X X 

Q = h = <p|z|q>  , 

pq  pq 

since  the  AO  basis  functions  are  not  dependent  upon  an  electric  field 
perturbation.  Also,  for  an  electric  field  perturbation,  the  total 
derivative  given  in  Equation  [28]  reduces  to  the  molecular  orbital 
contribution;  the  AO  derivative  part  is  zero,  and  the  remaining  terms 
vanish  because  the  elements  of  are  zero.  The  perturbation  may  be 

the  delta  function. 


e 

X 0 = X s S(r.  - r)  . 
i 

If  so,  the  energy  derivative  is  the  electron  density  evaluated  at 
point  r.  The  general  perturbation  matrix  for  a delta  function  is 
given  by: 
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X X 

Q - h » <p|  6(r'-r)  |q>  = • 

pq  pq 


At  the  Hartree-Fock  level  of  theory,  response  properties  are 
equivalent  to  the  expectation  value  since  the  Hellman-Feynmann  theorem 
is  satisfied.  This  is  commonly  expressed  as 


SCF 

9E  SCF  X 

-EDO 

ax  X=0  y,v  yv  yv  , 

where  is  the  SCF  density  matrix  in  the  primitive  atomic  orbital 

basis. 

Analytic  response  properties  corresponding  to  the  derivative  of 
the  correlation  energy  are  evaluated  in  an  expression  similar  to 

the  SCF  energy  derivative, 


total 

3E  total  X 

-ED  Q 

ax  X=0  y,v  yv  yv  , 

where  the  SCF  density  is  corrected  by  correlation  contributions 

total  SCF 

D = D + D 

yv  yv  yv  . 

Returning  to  Equation  (15),  the  derivative  of  the  MBPT/CC  energy  is 
expressed  as  the  sum  of  atomic  and  molecular  orbital  contributions. 

If  the  basis  functions  remain  stationary  and  are  independent  of  the 
perturbation,  the  atomic  orbital  contribution  is  zero.  The  remaining 
contributions  may  be  expressed  as  follows: 
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a,b 


+ 2 £ D .Q-^. 
. ai  ai 
a,i 


or 


AE^ 


p.q 


D 

pq  pq 


In  addition,  for  the  delta-function  perturbation  the  derivative  energy 
becomes: 


^ ■’pq  ♦q'"*  ’ 

p,q 

Therefore  the  matrix  D may  be  identified  as  the  matrix  representation  of 
the  correlation  correction  to  the  electron  density.  The  equations  for  the 
density  matrix  may  be  easily  identified  from  Tables  (8-13)  where  is 
identically  zero. 

The  solution  of  the  CPHF  equation  introduces  the  full  effect  of  orbital 
relaxation  resulting  in  an  infinite-order  summation  of  certain  double- 
perturbation diagrams  [48]  and  should  be  preferable  to  any  average  value 
approach  based  on  this  fact  alone.  In  addition,  the  response  density  is  a 
by-product  of  an  analytic  gradient  calculation.  The  computational 
convenience  of  such  approach  has  also  been  recognized  by  Simandiras,  et  al. 
in  their  route  to  MBPT(2)  dipole  derivatives  [51]. 


IMPLEMENTATION  OF  MANY-BODY  FIRST  DERIVATIVES 


The  final  expressions  as  defined  by  Equations  (31)  are  in 
computational  form.  The  implementation  of  gradients  for  CCD  and  MBPT 
models  require  the  following  steps: 

(a)  Solve  the  unperturbed  problem. 

(b)  Solve  the  A equation  if  necessary  (i.e.  this  is  necessary  for 
coupled-cluster  derivatives) 

(c)  Obtain  S^  and  for  each  degree  of  freedom. 

(d)  Construct  the  T intermediates  and  back-transform  to  the  primitive 
basis  set. 

(e)  Obtain  the  derivative  integrals  and  perform  a dot-product  with 
T(vvpa).  The  result  is  the  atomic  orbital  contribution. 

(f)  Construct  the  occupied-occupied  and  virtual-virtual  blocks  of  the 
density. 

(g)  Form  the  occupied-occupied  and  virtual-virtual  blocks  of  the  I 
intermediate  utilizing  T. 

(h)  Solution  of  the  X and  I virtual-occupied  contributions. 

(i)  Solve  the  linear  equation  given  in  Equation  (30)  to  form  the 
virtual-occupied  block  of  the  density  matrix. 

(j)  Perform  necessary  dot  products  for  each  3N  degrees  of  freedom  and 
add  the  corresponding  atomic  orbital  contribution. 

(k)  Compute  all  desired  one-electron  properties. 

Although  the  density  matrix  is  constructed  as  a by-product  of  a 
gradient  calculation,  it  is  sometimes  desirable  to  calculate  it 
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independently.  For  a given  MBPT/CC  model,  the  occupied-occupied  and 
virtual-virtual  blocks  of  D and  the  intermediate  X are  constructed. 
Equation  (30)  is  then  solved  for  the  virtual-occupied  block  of  D.  The 
SCF  density  in  the  molecular  orbital  basis  “^ij^  added  to 

obtain  the  total  density  matrix  in  the  molecular  orbital  basis  set. 

The  density  matrix  is  then  transformed  to  the  atomic  orbital  basis  so 
that  response  properties  may  be  evaluated  as  in  Equation  (31).  Thus 
the  transformation  of  each  set  of  one-electron  property  integrals  to 
the  molecular  orbital  basis  is  avoided;  only  the  transformation  of  the 
density  is  necessary.  The  construction  of  the  occupied-occupied  an 
virtual-virtual  blocks  of  the  density  matrix  and  the  intermediate  X 
are  for  the  MBPT(2),  MBPT(3),  MBPT(4)  and  CCD  models  are  identified  in 
Tables  (8-13).  Each  set  of  equations  is  complete  through  the  given 
order  in  correlation. 

There  are  several  relationships  that  may  be  exploited  in  the  route 
to  efficient  many-body  first  derivatives.  Naturally,  the  most 
Important  savings  are  due  to  the  ability  to  solve  perturbation 
independent  equations  thus  avoiding  derivative  t-amplitudes  and 
ultimately  leading  to  a response  density.  These  points  have  already 
been  discussed  at  some  length.  There  remains,  however,  several 
computational  simplifications  that  should  be  discussed. 

Finite-order  MBPT  derivatives  beyond  MBPT(2)  are  in  some  sense 
more  tedious  than  CC  expressions.  This  statement  is  easily  justified 
by  comparing  the  expressions  for  MBPT  first  derivative  given  in  Tables 
(10-12)  vs  the  CCD  derivative  expressions  in  Tables  (12-13). 

Naturally  the  inclusion  of  Tj^  and  T^  terms  complicate  the  MBPT(4) 
expression,  but  the  complexity  of  implementing  MBPT  derivatives 
arising  from  the  need  to  explicitly  consider  the  iterates  of  the  t- 
amplitudes. 
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An  excellent  example  is  the  construction  of  the  occupied-occupied 
and  virtual-virtual  blocks  of  the  density  matrix.  In  the  case  of  the 
CCD  model,  this  block  is  merely  the  product  of  converged  t2~  and 
X2-amplitudes.  By  contrast  the  occupied-occupied  block  for  MBPT(4) 
appears  to  require  six  such  products.  However,  if  one  manipulates  the 
expressions,  this  block  is  reduced  to  four  terms.  The  t^^  and  t^  terms 
are  evaluated  as  given  in  Table  (11),  but  the  terms  involving  are 
computed  by; 

ab  ab  ab  ab 

t (1)  t + At  (2)  At  (2) 
jk  ki  jk  ki 

where 

ab  ab  ab  ab  ab 

t = t (1)  + 2At  (3)  + 2At  (3)  + 1/2  Q (t,txt)  . 
ij  ij  ij  ij  ij 

The  term  t bar  is  easily  formed  from  r(ij,ab).  The  virtual-virtual 
block  is  also  constructed  using  t bar. 

The  form  of  the  diagonal  blocks  for  fourth-order  raises  another 
interesting  point.  Inspection  of  the  terms  that  contribute  to  the 
intermediate  X reveals  four  terms  that  are  exactly  the  diagonal  blocks 
contracted  with  an  integral;  however,  for  any  case  where  the  occupied- 
occupied  and  virtual-virtual  blocks  are  symmetric  only  two  of  the 
terms  need  to  be  calculated  (see  Table  (8)).  Naturally,  the 
introduction  of  A causes  these  blocks  to  be  inherently  unsymmetric, 
but  one  has  the  flexability  to  force  these  blocks  to  be  symmetric 
since  any  dot  product  performed  with  the  total  density  involves  a 
symmetric  matrix. 
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In  the  general  expressions,  the  products  that  form  the  I and  X 

intermediates  may  also  be  simplified  by  a similar  analysis.  Here,  The 

simplification  may  be  carried  out  if  the  terms  of  the  form  T(pq,rs) 

are  of  a proper  symmetry.  Consider  an  example,  the  contribution  to 

I . of  the  form: 
ai 

-1/2  E {r(ab,cd)  - r(ba,cd)  + T(cd,ab)  + r(cd,ba))  <ib||cd> 
b,c,d 

may  be  simplified  to 

-2  E r(ab,cd)  <ib| |cd> 
b,c,d 

provided  that  r(ab,cd)  is  antisymmetric  for  the  labels  a,b  and  c,d  and 
is  symmetric  to  interchange  of  the  pair  of  labels  ab  and  cd.  These 
conditions  are  met  for  this  term  and  all  other  such  terms  in  Table  (8) 
for  perturbation  theory.  This  follows  from  the  form  of  the  F 
intermediates  and  the  fact  that  the  t-amplitudes  are  antisymmetric. 
This  is  a fairly  significant  simplification  since  it  eliminates 
twenty-four  terms  where  N is  the  number  of  basis  functions.  The 
simplifications  due  to  symmetric  pairwise  interchanges  (interchanging 
labels  between  t-  and  X-amplitudes)  do  not  occur  for  CCD;  however,  the 
simplifications  may  be  made  by  imposing  symmetry  for  these  type  of 
interchanges  onto  the  F intermediates.  If  the  F intermediates  are 
averaged 

_ tu  tu  tu  tu 

F(pq,rs)  =1/2  E ( t X + X t ) 
t,u  pq  rs  pq  rs 

then  there  is  a savings  of  twelve  terms,  but,  if  done 
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straightforwardly,  at  the  expense  of  an  additional  step  for  each  T; 
this  quickly  becomes  prohibitive  when  one  considers  spin  cases  for 
each  r.  Fortunately,  this  is  not  the  only  way  to  symmetrize  F.  The 
symmetric  gamma  may  also  be  constructed  by  constructing  one  half  of  F 
as  it  appears  in  Table  (13),  sorting  the  labels  and  adding  the  result. 
This  is  at  a cost  of  only  N^.  Exactly  the  same  type  treatment  may  be 
used  for  the  construction  of  F terms  arising  from  products  of 
t2-amplitudes  for  MBPT(4). 

The  above  simplifications  also  have  the  effect  of  allowing  one  set 
of  routines  to  build  F and  to  perform  the  necessary  products.  Such 
relationships  are  inherent  in  the  comparison  of  MBPT  and  CC.  A well 
written  CC  code  should  exploit  the  relationship  with  finite-order 
MBPT. 

The  terms  that  arise  from  quadruple  excitations  in  the  derivative 
expressions  are  closely  related  to  the  normal  quadratic  energy  term. 

In  fact,  all  new  terms  arising  from  quadruples  may  be  determined  using 
existing  computational  procedures.  This  requires  the  following 
permutations  of  the  function  Q(I,txt);  (1)  Q(X,txt),  (2)  Q(t,IxX),  and 
(3)  Q(t,XxI).  Only  the  first  term  is  needed  for  fourth-order  as 
demonstrated  in  the  Appendix.  The  second  and  third  terms  are  used  in 
the  evaluation  of  A for  the  CCD  model. 

If  one  recalls  that  the  derivative  expressions  differ  only  by  the 
quadruple  excitation  terms,  there  is  justification  in  using  the 
converged  t-amplitudes  as  the  initial  approximation  to  A.  In  practice 
the  coupled-cluster  equations  are  solved  using  the  reduced  linear 
equation  method  (RLE)  to  accelerate  convergence  [52].  The  same 
procedure  is  employed  for  solving  the  linear  equation  A.  Studies 
show  that  CC  equations  converge  in  10-16  iterations  using  RLE  [52,53]. 
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The  experience  with  the  CCD  gradient  code  suggests  that  the  linear  A 
equation  converges  in  approximately  one-third  to  one-half  the  number 
of  iterations  required  to  solve  the  unperturbed  problem. 

It  should  be  noted  that  the  fourth-order  energy  gradient  is  in 
the  number  of  basis  functions.  This  dependence  is  due  to  the  triple- 
excitation terms,  which  are  an  procedure  in  the  unperturbed  energy 
expressions.  However,  in  addition  to  the  calculation  of  T^ 
amplitudes,  the  T^  contribution  to  T2,  the  contribution  to  the 
diagonal  blocks  of  the  density,  and  the  product  of  T2  and  T^  needed  to 
form  r(ai,bc)  and  r(jk,ia)  introduce  costly  steps.  The  computational 
overhead  is  close  to  that  of  determining  the  energy  contribution  from 
triples  in  a single  iteration  of  the  CCSDT-1  model.  The  reason  for 
this  lies  in  the  fact  that  the  fourth-order  energy  may  be  calculated 
by  the  formula 

T abc  2 

E = E (t  ) x(e+e+e-e-e-e), 

4 ijk  i j k a b c 

a>b>c 

i>j>k 

thereby  eliminating  the  need  for  the  inclusion  of  T^  to  T2«  This 
symmetry  does  not  hold  for  the  case  of  energy  derivatives.  If  one 
considers  that  a numerical  calculation  of  a single  property  typically 
requires  at  least  two  calculations,  the  additional  computational 
overhead  of  an  analytic  calculation  is  more  than  justified. 

It  should  be  pointed  out  that  all  computer  codes  were  thoroughly 
checked  by  comparisons  with  numerically  calculated  one-electron 
properties  and  gradients.  Finite-field  calculations  were  normally 
reproduced  to  six  decimal  places  which  represents  the  expected 
accuracy  for  numerical  results. 


APPLICATIONS  OF  MANY-BODY  FIRST  DERIVATIVES 


Stationary  Properties 

Bauschlicher  and  Taylor  published  a recent  report  featuring  full 
Cl  polarizabilities  and  dipole  moments  for  a series  of  small  systems 
[54].  Whereas  the  basis  sets  used  are  not  adequate  for  quantitative 
comparison  with  experiment,  such  calculations  do  provide  important 
benchmarks  for  the  calibration  of  the  various  treatments  of  electron 
correlation.  This  situation  is  not  unlike  full  Cl  benchmarks  for 
total  energies.  These  Cl  results  are  available  due  to  improvements  in 
methodology  and  computer  technology  (These  calculations  were  performed 
on  a CRAY  2 featuring  256  Megawords  of  main  memory). 

In  this  work,  CCD  and  SDTQ-MBPT(4)  polarizabilities  and  dipole 
moments  are  compared  with  those  obtained  by  full  Cl  and  other  levels 
of  theory  (54)  for  the  state  of  HF,  the  A^  state  of  CH2,  and  F . 
In  addition,  other  CCD  and  MBPT(4)  one-electron  properties  for  these 
systems  are  tabulated. 

The  [4s2pld]/[2slp]  basis  sets  for  HF  and  CH2,  the  [5s4p2d]  basis 
set  for  F“  and  the  molecular  geometries  are  those  reported  in 
Reference  [54].  Unlike  Bauschlicher  and  Taylor's  treatment  of  the 
basis  sets,  the  3s  combination  of  the  3d  functions  has  not  been 
eliminated.  In  addition,  the  wavefunction  has  been  improved  by 
correlating  all  electrons.  The  difference  due  to  the  3s  component  is 
in  the  millihartree  range  or  less  at  the  SCF  level  of  theory.  In  the 
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case  of  HF,  where  calculations  with  a frozen  core  and  the  3s 
combination  eliminated  are  available  [21],  the  MBPT(4)  energies  are 
lower  by  0.02  hartrees.  Nevertheless,  a comparison  of  dipole  moments 
and  polarizabilities  is  valid  since  the  core  region  is  not  expected  to 
make  significant  contributions  to  either  property. 

The  dipole  moment  and  polarizability  are  computed  along  the  bond 
axis  for  HF  and  the  axis  for  CH2.  All  MBPT(4)  and  CCD  first-order 
properties  are  computed  analytically.  The  polarizability  is  computed 
by  solving  for  the  response  density  in  the  presence  of  an  external 
electric  field. 

The  results  for  CH2  are  summarized  in  Table  (15).  The  state 
of  CH2  is  not  adequately  described  by  a single  reference.  The  MBPT(4) 
values  are  in  good  agreement  with  the  SDCI,  yet  both  are  a few  percent 
from  the  full  Cl  values.  Moreover,  the  CCD  values  are  slightly 
Improved  over  both  MBPT(4)  and  SDCI.  This  suggests  that  the  inclusion 
of  higher  excitations  are  necessary  to  overcome  the  poor  reference 
state.  Additional  properties  of  CH2  at  the  MBPT(4)  and  CCD  levels  of 
theory  are  listed  in  Table  (16).  The  majority  of  properties  are  only 
slightly  changed  with  the  addition  of  correlation  with  one  of  the 
exceptions  being  the  dipole  moment. 

The  polarizability  for  F“  is  given  in  Table  (17).  There  is  a 
large  change  for  any  method  relative  to  SCF.  The  MBPT(4)  value  is 
1.85  times  the  SCF  polarizability.  The  MBPT(4)  value  overshoots  the 
full  Cl  while  the  SDCI  and  CCD  values  are  below  by  a slightly  larger 
amount.  Cole  and  Bartlett's  [21]  comparison  of  the  MBPT(4)  energies 
with  full  Cl,  using  comparable  basis  sets,  shows  the  fourth-order 
triples  cause  the  MBPT  energy  to  overshoot  that  of  the  full  Cl.  They 
conclude  that  this  is  due  to  the  relative  importance  of  higher 


excitations  in  F . 
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Three  bond  lengths  were  considered  for  HF.  Comparisons  are  given 
in  Tables  (18-20).  The  SCF  reference  falls  as  the  bond  is  stretched 
and  correlation  becomes  more  important;  this  is  reflected  in  the  list 
of  properties  for  HF  given  in  Tables  (21-23).  At  equilibrium,  the 
agreement  of  CCD  and  MBPT(A)  with  the  full  Cl  is  excellent.  The 
MBPT(4)  polarizability  is  the  same  as  the  SDCI+Q  (SDCI  with  the 
inclusion  of  the  Davidson  correction)  result,  and  the  dipole  moment  is 
close  to  the  multireference  values.  The  CCD  results  are  extremely 
good;  the  polarizability  is  equal  to  that  obtained  by  the  full  Cl  and 
the  dipole  moment  differs  by  less  than  one  percent.  At  1.5  x r^,  the 
CCD  and  MBPT(4)  dipole  moments  are  closer  to  the  full  Cl  that  that 
obtained  by  SDCI.  The  fourth-order  polarizability  is  lower  than  all 
other  methods,  but  it  is  closer  to  the  full  Cl  than  the  SDCI+Q  value. 
The  CCD  and  MBPT(4)  dipole  moments  are  substantially  better  than  that 
obtained  by  the  CISD  or  CISD+Q  at  2 x r^.  The  MBPT(4)  polarizability 
is  equal  to  that  obtained  by  the  SDCI+Q;  however,  the  polarizability 
obtained  by  CCD  is  surprisingly  high  suggesting  that  the  additional 
correlation  effects  due  to  single  excitations  are  important. 

In  conclusion,  the  overall  comparison  of  CCD  and  SDTQ-MBPT(4) 
properties  to  those  obtained  by  full  Cl  is  very  favorable.  It  is 
interesting  to  note  how  well  the  fourth-order  CCD  and  MBPT  results 
parallel  the  truncated  Cl  for  the  cases  presented.  Comparison  of  CCD 
and  MBPT(4)  suggests  that  infinite  order  models  are  superior  for 
second  order  properties  such  as  the  polarizability. 

An  in-depth  property  study  of  H2O  using  the  84-function  basis  set 
of  Davidson  and  Feller  [55]  has  been  performed.  A comparison  of  MBPT 
and  CC  response  properties  with  SCF,  Cl  and  experimental  properties  is 
given  in  Table  (24).  Often  a large  basis  set  is  needed  to  adequately 
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describe  one-electron  properties.  When  a specific  property  Is 
desired,  the  basis  set  may  be  chosen  to  describe  the  important  regions 
associated  with  the  property.  This  basis  set  provides  excellent 
agreement  with  the  full  spectrum  of  properties  considered.  Feller,  et 
al.  [56]  performed  a study  on  the  convergence  of  this  basis  set  and 
found  that  there  were  only  minor  changes  due  to  increasing  the  size  of 
the  basis  set. 

It  should  be  pointed  out  that  the  response  density  method  allows 
for  such  a large  scale  property  study;  if  traditional  finite-field 
techniques  were  employed  this  study  would  be  nearly  prohibitive  due  to 
the  necessary  number  of  calculations.  In  order  to  compute  the 
properties  given  in  Table  (24),  31  finite-field  calculations  would  be 
necessary  for  each  level  of  theory.  Analytic  techniques  reduce  the 
number  of  calculations  to  one  per  level  of  theory. 

In  the  past,  attempts  to  compute  analytical  MBPT/CC  properties 
have  been  limited  to  expectation  value  approaches  [50,57,58].  Table 
(24)  also  contains  MBPT(2)  properties  computed  as  an  expectation 
values  of  the  wavefunction.  The  response  properties  have  the 
advantage  of  summing  many  higher-order  diagrams.  Consequently,  the 
energy  derivative  expression  should  be  preferable  to  the  expectation 
value  approach  on  that  basis  alone. 

The  general  trends  of  the  effects  of  correlation  Davidson  and 
Feller  reported  are  not  changed  here;  however,  the  inexpensive  MBPT 
calculations  produce  improvement  over  the  Cl  wavefunction  for  some 
properties,  specifically  the  dipole  moment  and  oxygen  field  gradient. 
It  is  interesting  to  note  that  the  many-body  results  are  closer  in 
agreement  with  the  multireference  Cl  that  the  single  reference  Cl. 

The  overall  comparison  with  experiment  is  ±2-3  percent. 
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If  the  method  exactly  satisfied  the  Hellman-Feynmann  theorem,  and 
if  the  molecule  were  at  the  predicted  equilibrium  geometry,  the  forces 
would  vanish.  Hence,  the  slightly  larger  value  here  can  be  attributed 
to  either  the  Cl  geometries  being  closer  to  the  experimental  value  or 
the  non-satisfaction  of  the  Hellman-Feynmann  theorem. 

A measure  of  the  orbital  relaxation  in  first-order  properties  is 
given  by  the  difference  between  the  results  in  MBPT(2)^^g  (Equation 
(3A))  and  the  corresponding  response  density,  which  includes  the 
relaxation  effects  to  all  orders,  while  the  correlation  effect  is  only 
second  order.  The  average  effect  is  about  1.5  percent.  A MBPT(4)g^^g 
density  will  have  recovered  more  relaxation  terms  making  it  nearly 
equivalent  to  the  corresponding  MBPT(4)  response  density.  The  reason 
relaxation  effects  are  relatively  small  here  is  that  the  first-order 
property  is  based  upon  the  response  of  the  SCF  wavefunction, 
therefore,  orbital  relaxation  is  only  introduced  via  coupling  through 
double-excitation  correlation  corrections.  However,  in  a second-order 
property  like  a polarizability,  one  might  expect  significant 
differences  due  to  orbital  relaxation  [59,60]. 

Molecular  Gradients 


The  development  of  molecular  gradients  created  a viable  means  to 
calculate  both  theoretical  force  constants  and  molecular  geometries. 
The  force  method  developed  by  Pulay  [2]  for  use  at  the  SCF  level  of 
theory  suggested  that  force  constants  be  calculated  by  first 
analytically  differentiating  the  energy  followed  by  numerical 
differentiation.  The  method  has  proven  to  be  reliable  and  has  been 
extended  to  correlated  approaches  such  as  Cl  [3],  MCSCF  [4],  MBPT(2) 
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[5],  and  CC  122,23,36].  In  this  work  this  method  is  applied  to  SDTQ- 
MBPT(4)  for  the  first  time.  In  addition,  initial  applications  of  an 
efficient  procedure  for  CCD  is  presented. 

The  applications  chosen  here  were  based  on  the  recent  reports  of 
CCSD  and  CCSDT-1  harmonic  vibrational  frequencies  [22,23,60]. 
Configuration  Interaction  results  are  also  available  [60,61].  The 
initial  applications  here  are  meant  to  test  the  reliability  of  the 
SDTQ-MBPT(4)  and  CCD  models  in  comparison  with  Cl,  CCSD  and 
experiment. 

In  order  to  provide  direct  comparisons,  the  basis  sets  employed  in 
this  study  are  those  used  in  Reference  [60].  These  are  standard 
double  zeta  plus  polarization  (DZP)  Huzinaga-Dunning  basis  sets 
designated  as  (9s5pld/4s2pld)  [62,63].  Polarization  function  exponents 
were  a^(C)  - 0.75,  a^(N)  - 0.80,  a^(0)  = 0.85  and  ap(H)  = 1.0.  The 
hydrogen  s functions  were  scaled  by  a factor  of  1.2. 

Tables  (25,26)  contain  molecular  constants  calculated  by  both  CCD 
and  MBPT(4)  for  H2O  and  NH^.  Tables  (27,28)  contain  CCD  molecular 
constants  for  HCN  and  CH^.  Included  are  equilibrium  geometries, 
harmonic  vibrational  frequencies,  and  dipole  moments  for  several 
levels  of  theory. 

H^O.  The  geometry  of  H2O  is  well  described  by  all  levels  of 
theory  as  compared  with  experiment.  Harmonic  vibrational  frequencies 
steadily  improve  as  one  follows  the  progression  of  CCD,  CCSD,  MBPT(4) 
and  CCSDT-1  except  for  the  «2  where  the  MBPT(4)  frequency  is  slightly 
lower  than  that  obtained  by  CCSDT-1.  The  improvement  due  to  triple 
excitations  is  consistent  with  a study  by  Bartlett,  et  al.  of  the 
quartic  force  field  of  H2O.  A comparison  with  the  CISDTQ  method  is 
also  valid  since  it  is  expected  to  be  close  to  the  full  Cl  for  a ten 
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electron  system.  The  discrepancies  of  CCD  and  MBPT(4)  with  CISDTQ  are 
36,  8,  35  and  14,  3,  14,  respectively.  Overall  the  average  error  of 
MBPT(4)  frequencies  compared  with  experiment  is  approximately  two 
percent  while  those  for  CCD  are  around  2.5  percent.  The  dipole  moment 
follows  the  same  basic  trend  as  the  reported  frequencies;  however, 
this  relatively  small  basis  set  is  inadequate  for  a valid  comparison 
with  experiment;  this  should  be  obvious  from  the  earlier  discussion  of 
H2O  properties. 

NH^.  The  geometry  is  essentially  the  same  from  all  correlated 
treatments  and  is  in  good  agreement  with  experiment.  The  totally 
symmetric  mode  of  ammonia  is  in  error  for  all  reported  theoretical 
techniques.  The  best  theoretical  value  is  the  MBPT(4)  result  with  a 
average  error  compared  with  experiment  of  8.51  percent.  It  has  been 
suggested  [60]  that  this  problem  is  related  to  the  small  basis  set 
used  in  this  study.  A comparison  with  the  CISDTQ  result  reveals  that 
MBPT(4),  CCD,  CCSD  and  CCSDT-1  approach  the  basis  set  limit. 

HCN.  The  geometry  and  dipole  moment  of  HCN  is  best  described  by 
the  CISD  followed  by  CCD  then  CCSD.  This  trend  is  exactly  reversed 
for  harmonic  frequencies.  The  average  error  is  4.16  percent,  1.51 
percent  and  0.82  percent  for  CISD,  CCD  and  CCSD,  respectively. 

CH,.  The  results  for  methane  parallel  the  result  for  HCN. 

— 4 

Coupled-cluster  methods  are  closer  to  experiment  than  CISD,  but  the 
differences  are  not  as  exaggerated  as  for  HCN.  The  differences 
between  CCD  and  CCSD  are  7 cm"^  (w^),  2 cm"^  ((O2),  11  cm"^  (w^),  and  2 
cm“^  respectively.  The  average  error  is  2.37  percent,  1.67 

percent,  and  1.45  percent  for  CISD,  CCD  and  CCSD,  respectively. 

These  initial  investigations  suggest  that  many-body  methods  are 
well-suited  for  the  description  of  harmonic  vibrational  frequencies 
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and  molecular  structures.  In  all  cases,  the  MBPT/CC  results  are 
significantly  improved  over  truncated  Cl.  The  comparison  of  CCD  with 
CCSD  suggests  that  CC  theory  rapidly  converges  for  the  systems 
studied.  It  is  also  encouraging  to  find  that  the  non-iterative 
MBPT(4)  results  closely  parallel  CCSDT-1.  If  this  trend  continues 
upon  further  investigation,  MBPT(4)  will  provide  an  important  route 
for  geometry  optimizations  and  prediction  of  harmonic  frequencies  for 
extended  basis  studies. 
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Table  15.  Properties 

for  CH2 

Level  of  Theory 

Energy 

a 

y 

SCF(5d) 

-38.886296 

10.2 

0.0807 

SCF  (6d) 

-38.886297 

10.1 

0.807 

MBPT(4) 

-39.036902 

10.4 

0.754 

CCD 

-39.037433 

10.5 

0.734 

FCI 

-39.027182 

10.7 

0.716 

SDCI 

-39.018283 

10.4 

0.743 

SDCI+Q 

-39.027221 

10.7 

0.743 

CPF 

-39.025116 

10.7 

0.721 

2CPF 

-38.907659 

10.6 

0.716 

2CPF+CI 

-39.022155 

10.7 

0.717 

2CPF+CI+Q 

-39.027741 

10.7 

0.715 

CASSCF 

-38.945528 

10.4 

0.691 

MRCI 

-39.025803 

10.7 

0.713 

MRCI+Q 

-39.028479 

10.7 

0.716 

Basis  set  and  geometry  are  defined  in  Reference  [54]. 
CI,CASSCF,  CPF  and  multireference  results  taken  from 
Reference  [54].  All  dipole  moments  computed  as  energy 
derivatives.  All  properties  are  given  in  atomic  units. 
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Table  16.  Comparison  of  SCF  and  Many-Body  Properties  for  CH2 


Property 

SCF 

SDTQ-MBPT(4) 

CCD 

Energy 

-38.886297 

-39.036902 

-39.037433 

0.8070 

0.7543 

0.7348 

0(xx) 

0.3970 

0.3928 

0.3858 

e(yy) 

1.2285 

1.1807 

1.1249 

0(zz) 

-1.6255 

-1.5735 

-1.5107 

-15.6399 

-15.6463 

-15.6474 

(l/rH)e- 

-4.2320 

-4.2340 

-4.2365 

qz/rl 

-0.0685 

-0.0670 

-0.0648 

3 

qz/rjj 

0.0032 

0.0091 

0.0084 

qy/r„ 

0.0017 

0.0142 

0.0141 

qxxc 

-0.4304 

-0.4422 

-0.4329 

qyyc 

1.4984 

1.4231 

1.3810 

qzzc 

1.0680 

0.9810 

0.9489 

qxxjj 

0.1077 

0.1108 

0.1113 

qyyfl 

-0.1005 

-0.1052 

-0.1052 

qzzjj 

-0.0072 

-0.0057 

-0.0052 

All  properties  given  in  atomic  units.  Properties  are  calculated 
relative  to  the  origin  unless  otherwise  denoted. 
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Table  17.  Properties  for  F 


Level  of  Theory 

Energy 

a 

SCF(5d) 

-99.444462 

9.894 

SCF  (6d) 

-99.446023 

9.909 

MBPT(4) 

-99.714359 

18.362 

CCD 

-99.697915 

13.157 

FCI 

-99.659493 

16.295 

CASSCF(2p) 

-99.544014 

16.903 

CASSCF(2s2p) 

-99.582645 

13.906 

SDCI 

-99.645888 

13.965 

SDCI+Q 

-99.656183 

15.540 

CPF 

-99.654584 

16.050 

MRCI 

-99.657746 

16.134 

MRCI+Q 

-99.660278 

16.346 

MRCIBIG 

-99.658584 

16.034 

MRCIBIG+Q 

-99.659593 

16.303 

Basis  set  and  geometry  are  defined  in  Reference  [54]. 
CI,CASSCF,  CPF  and  multireference  results  taken  from 
Reference  [54].  MRCIBIG  contains  single  and  double 
excitations  relative  to  the  CASSCF(2s2p)  reference. 

All  dipole  moments  computed  as  energy 

derivatives.  All  properties  are  given  in  atomic  units. 
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Table  18.  Properties 

for  HF 

Level  of  Theory 

Energy 

a 

V 

SCF(5d) 

-100.047086 

4.2 

0.814 

SCF(6d) 

-100.047620 

4.3 

0.812 

MBPT(4) 

-100.268771 

4.5 

0.762 

CCD 

-100.264540 

4.4 

0.769 

FCI 

-100.250968 

4.4 

0.765 

SDCI 

-100.241588 

4.4 

0.773 

SDCI+Q 

-100.249372 

4.5 

0.766 

CPF 

-100.247719 

4.4 

0.767 

CASSCF 

-100.070188 

4.4 

0.766 

MRCI 

-100.244725 

4.4 

0.765 

MRCI+Q 

-100.250863 

4.4 

0.763 

Basis  set  and  geometry  are  defined  in  Reference  [54]. 
CI,CASSCF,  CPF  and  multireference  results  taken  from 
Reference  [54].  All  dipole  moments  computed  as  energy 
derivatives.  All  properties  are  given  in  atomic  units. 
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Table  19.  Properties 

for  HF  1.5xr^ 

Level  of  Theory 

Energy 

a 

V 

SCF(5d) 

-99.933228 

12.2 

1.161 

SCF(6d) 

-99.933801 

12.2 

1.162 

MBPT(4) 

-100.176935 

11.9 

0.925 

CCD 

-100.169302 

12.6 

0.929 

FCI 

-100.160392 

12.2 

0.896 

SDCI 

-100.145522 

12.3 

0.951 

SDCI+Q 

-100.158979 

12.6 

0.879 

CPF 

-100.156162 

12.2 

0.900 

CASSCF 

-99.9864238 

12.4 

0.784 

MRCI 

-100.154120 

12.3 

0.890 

MRCI+Q 

-100.160276 

12.2 

0.896 

Basis  set  and  geometry  are  defined  in  Reference  [54]. 
CI,CASSCF,  CPF  and  multireference  results  taken  from 
Reference  [54].  All  dipole  moments  computed  as  energy 
derivatives.  All  properties  are  given  in  atomic  units. 
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Table  20.  Properties 

for  HF  2xr 

e 

Level  of  Theory 

Energy 

a 

y 

SCF(5d) 

-99.817571 

26.2 

1.511 

SCF(6d) 

-99.818241 

26.2 

1.516 

MBPT(4) 

-100.093551 

21.6 

0.717 

CCD 

-100.076194 

28.6 

0.652 

FCI 

-100.081107 

19.1 

0.618 

SDCI 

-100.053534 

22.8 

0.811 

SDCI+Q 

-100.082399 

21.6 

0.451 

CPF 

-100.075803 

19.7 

0.645 

CASSCF 

-99.920994 

13.7 

0.418 

MRCI 

-100.075566 

18.6 

0.598 

MRCI+Q 

-100.080738 

19.0 

0.615 

Basis  set  and  geometry  are  defined  in  Reference  [54]. 
CI,CASSCF,  CPF  and  multireference  results  taken  from 
Reference  [54].  All  dipole  moments  computed  as  energy 
derivatives.  All  properties  are  given  in  atomic  units. 
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Table  21.  Comparison  of  SCF  and  Many-Body  Properties  HF  (r^) 


Property 

SCF 

MBPT(4) 

CCD 

Energy 

-100.047620 

-100.268771 

-100.264540 

0.8116 

0.7617 

0.7687 

e(xx) 

-0.9112 

-0.8741 

-0.8811 

0(zz) 

1.8224 

1.7483 

1.7623 

(l/rF)e- 

-27.1746 

-27.1458 

-27.1432 

(l/rH)e- 

-6.1105 

-6.1175 

-6.1172 

qz/rj 

-0.0837 

-0.0849 

-0.0850 

3 

qz/r^ 

0.0101 

0.0425 

0.0398 

qxxp 

-1.4907 

-1.3859 

-1.3811 

qzzp 

2.9814 

2.7719 

2.7622 

qxXp 

-0.2842 

-0.2967 

-0.2966 

qzzj^ 

0.5685 

0.5934 

0.5933 

All  properties  given  in  atomic  units.  Properties  are  calculated 
relative  to  the  origin  unless  otherwise  denoted. 
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Table  22.  Comparison  of  SCF  and  Many-Body  Properties  HF  (l.Sxr^) 


Property 

SCF 

MBPT(4) 

CCD 

Energy 

-99.933801 

-100.176935 

-100.169302 

1.1616 

0.9253 

0.9289 

0(XX) 

-1.8904 

-1.6033 

-1.6106 

6(zz) 

3.7809 

3.2067 

3.2213 

(l/rF>e- 

-26.9906 

-26.9530 

-26.9466 

(l/r,)e- 

-4.2401 

-4.3044 

-4.3043 

qz/rp 

-0.0527 

-0.0484 

-0.0498 

qz/tj 

-0.1465 

-0.1129 

-0.1191 

qxxp 

-1.9960 

-1.9108 

-1.9625 

qzzp 

3.9920 

3.8215 

3.9250 

qxxjj 

0.0075 

-0.0052 

-0.0046 

qzzjj 

0.0149 

0.0103 

-0.0092 

All  properties  given  in  atomic  units.  Properties  are  calculated 
relative  to  the  origin  unless  otherwise  denoted. 
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Table  23.  Comparison  of  SCF  and  Many-Body  Properties  HF  (2xr^) 


Property 

SCF 

MBPT(4) 

CCD 

Energy 

-99.818241 

-100.093551 

-100.076194 

1.5165 

0.7168 

0.6524 

0(xx) 

-3.1218 

1.7493 

1.6768 

0(zz) 

6.2436 

3.4985 

3.3536 

<l/^F>e- 

-26.9113 

-26.8343 

-26.8096 

(l/r,)e- 

-3.2706 

-3.4568 

-3.4577 

, 3 

qz/rp 

-0.0335 

-0.0208 

-0.0255 

, 3 

qz/rjj 

-0.1043 

-0.0593 

-0.0761 

qxxp 

-2.1979 

-2.2945 

-2.5662 

qzzp 

4.3958 

4.5891 

5.1325 

qxxy 

0.0209 

0.0074 

0.0097 

qzZjj 

-0.0417 

-0.0149 

-0.0194 

All  properties  given  in  atomic  units.  Properties  are  calculated 
relative  to  the  origin  unless  otherwise  denoted. 
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Table  24. 

Properties  for  H_0 
Experimental  Geometry; 
in  atomic  units. 

84  GTOs;  all 

properties 

are  given 

Property 

SCF 

MBPT(2)^^g 

MBPT(2) 

MBPT(3) 

D-MBPT(4) 

Energy 

-76.06661 

-76.37699 

-76.37699 

-76.38126 

-76.38604 

0.7951 

0.7465 

0.7518 

0.7571 

0.7521 

e(yy) 

1.8772 

1.8709 

1.9001 

1.8674 

1.8661 

n(e) 

0.8805 

0.8685 

0.8715 

0.8772 

0.8757 

23.440 

23.405 

23.447 

23.445 

23.445 

5.7766 

5.7686 

5.7633 

5.7736 

5.7730 

qz/ro 

0.0015 

0.0000 

0.0058 

0.0048 

0.0053 

qz/r^ 

0.0126 

0.0033 

0.0033 

0.0001 

0.0018 

, 3 
qy/r„ 

0.0125 

0.0042 

0.0048 

0.0012 

0.0032 

'lo 

-1.8512 

-1.7459 

-1.6544 

-1.7346 

-1.7248 

n(qo) 

0.7952 

0.7824 

0.7699 

0.7754 

0.7722 

% 

0.4702 

0.4822 

0.4805 

0.4805 

0.4819 

n(qo) 

0.1289 

0.1229 

0.1228 

0.1242 

0.1241 

Basis  set  and  geometry  from  Reference  [55].  Cl  and  experimental 
values  taken  from  Reference  [55].  All  properties  are  defined  in 
Reference  [55]. 
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SDQ-MBPT(4) 

SDTQ-MBPT(4) 

CCD 

SD-CI 

MRSD-CI 

EXPT. 

-76.38367 

-76.39306 

-76.38191 

-76.36784 

-76.37362 

-76.4376 

0.7514 

0.7422 

0.7561 

0.7611 

0.7588 

0.7296 

1.8648 

1.8708 

1.8689 

1.8701 

1.8669 

1.96 

0.8732 

0.8711 

0.8760 

0.8744 

0.8768 

0.90 

23.446 

23.443 

23.446 

23.451 

23.449 

23.36 

5.7725 

5.7697 

5.7732 

5.7766 

5.7738 

5.769 

0.0053 

0.0059 

0.0051 

0.0041 

0.0029 

0.0 

0.0018 

0.0042 

0.0010 

0.0007 

-0.0002 

0.0 

0.0032 

-0.0059 

0.0021 

-0.0003 

-0.0012 

0.0 

-1.7272 

-1.7044 

-1.7260 

-1.7484 

-1.7488 

-1.665 

0.7747 

0.7692 

0.7764 

0.7823 

0.7784 

0.75 

0.4813 

0.4830 

0.4776 

0.4791 

0.4802 

0.4583 

0.1245 

0.1237 

0.1215 

0.1252 

0.1261 

0.1350 
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Table  25.  Comparison  of  molecular  constants  for  H2O. 


SCF® 

CISD® 

CISDT® 

CISDQ^ 

CISDTQ® 

Energy 

-76.04695 

-76.257507 

-76.260391 

-76.267091 

-76.269799 

■■oh  <*> 

0.9441 

0.9585 

0.9592 

0.9621 

0.9629 

HOH  angle 

106.61 

104.86 

104.74 

104.59 

104.45 

Harmonic  ^ 
freq.  (cm  ) 

4164 

3967 

3953 

3908 

3892 

«2(ai) 

1751 

1693 

1689 

1681 

1677 

«3(b2> 

4287 

4090 

4077 

4037 

4022 

y (D) 

2.180 

2.144 

2.136 

2.138 

2.130 

All  calculations  performed  at  the  theoretical  equilibrium  geometry  using  a 
DZP  basis  set. 

^Reference  [23]. 

Reference  [64]. 

^Reference  [65]. 
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MBPT(4) 

CCD 

CCSD® 

CCSDT-1® 

Expt . 

-76.267764 

-76.265580 

-76.266438 

-76.269661 

0.9622 

0.9611 

0.9618 

0.9626 

0.957^^ 

104.51 

104.69 

104.61 

104.46 

104.5'’ 

3906 

3928 

3912 

3896 

3832*’ 

1674 

1685 

1683 

1676 

I649'’ 

4036 

4057 

4041 

4026 

3943*’ 

2.130 

2.143 

2.139 

2.130 

1.855'^ 
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Table  26.  Comparison  of  molecular  constants  for  NH^ 


SCF® 

CISD® 

CISDT® 

CISDTQ® 

Energy 

-56.209680 

-56.402670 

-56.406165 

-56.415178 

(A) 

1.001 

1.013 

1.014 

1.017 

HNH  angle 

108.2 

106.7 

106.5 

106.3 

Harmonic  ^ 
freq.  (cm”''’) 

3724 

3589 

3575 

3528 

(02(ai> 

1114 

1120 

1122 

1121 

«3(e) 

3872 

3735 

3721 

3676 

“4(e) 

1302 

1727 

1720 

1706 

y (D) 

1.805 

1.853 

1.856 

1.862 

All  calculations  performed  at  the  theoretical  equilibrium 
geometry  using  a DZP  basis  set. 

?Reference  [23]. 

Reference  [66]. 

^Reference  [67]. 

"^Reference  [68]. 
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MBPT(4) 

CCD 

CCSD® 

CCSDT-1®  Expt. 

-56.429773 

-56.426250 

-56.426818 

-56.430550 

1.015 

1.015 

1.015 

1.016 

1.012^ 

106.5 

106.5 

106.5 

106.5 

106. 7^’ 

3556 

3562 

3551 

3533 

3506^^ 

1109 

1120 

1119 

1121 

1022'^ 

3710 

3713 

3700 

3683 

3577*^ 

1710 

1718 

1716 

1709 

1691^^ 

1.851 

1.858 

1.850 

1.862 

1.472^^ 
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Table  27.  Comparison  of  molecular  constants  for  HCN. 


SCF® 

CISD® 

CCD 

CCSD® 

Expt.*^ 

Energy 

-92.88989 

-93.18841 

-93.20971 

-93.21288 

<*) 

1.0622 

1.0666 

1.0698 

1.0704 

1.065 

1.1368 

1.1618 

1.1686 

1.1709 

1.153 

Harmonic  ^ 

freq.  (cm  ) 

Wj^(a) 

3636 

3549 

3502 

3493 

3440 

858 

762 

734 

726 

727 

(03(a) 

2404 

2224 

2166 

2145 

2128 

y (D) 

3.219 

3.014 

2.956 

2.953 

2.985 

All  calculations  performed  at  the  theoretical  equilibrium  geometry  using  a 
DZP  basis  set. 
fpeference  [60]. 

Reference  [69]. 
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Table  28.  Comparison  of  molecular  constants  for  CH^. 


SCF® 

CISD^ 

CCD 

CCSD® 

Expt.*^ 

Energy 

-40.20759 

-40.39741 

-40.40548 

-40.40599 

■'CH 

1.0848 

1.0886 

1.0905 

1.0909 

1.0858 

Harmonic  ^ 

freq.  (cm  ) 

«j^(a) 

3174 

3116 

3094 

3087 

3026 

«2(e) 

1671 

1601 

1588 

1586 

1583 

»3(t2) 

3294 

3257 

3241 

3230 

3157 

1457 

1397 

1387 

1385 

1367 

All  calculations  performed  at  the  theoretical  equilibrium  geometry  using  a 
DZP  basis  set. 

?Reference  [60]. 

^Reference  [69]. 


CONCLUSIONS 


The  theory  and  implementation  of  CCD  and  MBPT  first  derivatives 
has  been  presented.  These  methods  allow  the  efficient  computation  of 
a variety  of  atomic  and  molecular  properties.  Initial  applications, 
intended  for  the  calibration  of  these  methods,  have  proven  to  be  in 
excellent  agreement  with  experiment  and  other  reliable  theoretical 
methods.  In  addition  to  the  work  presented  here,  the  graphical 
representation  of  the  electron  density  has  been  used  to  study  the 
dissociation  of  methane  [72]. 

The  methods  presented  here  are  currently  being  applied  to  a 
variety  of  systems.  Besler  et.  al  [60]  have  extended  their  studies  of 
CCSD  frequencies  to  nine  systems.  They  point  out  that  these  represent 
"a  majority  of  the  polyatomic  molecules  for  which  harmonic  vibration 
frequencies  are  known."  Many-body  studies  including  the  effects  of 
triple  excitations  are  limited  to  those  examples  discussed  here; 
however,  this  method  is  also  being  applied  to  other  systems.  Harmonic 
vibrational  frequencies  allow  for  a direct  comparison  of  theory  and 
experiment.  However,  these  quantities  are  difficult  to  deduce  from 
experiment.  Nevertheless,  the  information  obtained  in  these  studies 
will  allow  for  the  treatment  of  more  sophisticated  studies  including 
the  effects  of  anharmonicity.  A study  of  the  enharmonic  vibrational 
frequencies  for  methane  is  underway. 


103 


104 


The  theoretical  treatment  of  Infrared  Intensities  is  also  readily 
accessible  using  these  methods.  The  theoretical  study  of  intensities 
is  rarely  in  good  agreement  with  experiment.  The  reason  for  this 
discrepancy  is  attributed  to  the  uncertainty  of  experiment,  the  double 
harmonic  approximation  and  the  sensitivity  of  Intensities  to  the 
quality  of  the  basis  set. 

The  theory  present  here  has  been  presented  for  CCSDT  and  the  newly 
introduced  Quadratic  Configuration  Interaction  (QCI)  (73].  In 
addition,  the  techniques  and  algorithms  developed  for  MBPT/CC  are 
being  applied  to  a coupled-cluster  expectation  value  approach  termed 
XCC  (74). 

The  theory  of  analytic  second  derivatives  has  been  developed 
through  CCSDT-1  (73] . Analytic  second  derivatives  have  been  developed 
and  implemented  for  MBPT(2)  [75,76].  Efforts  to  Implement  higher 
levels  of  theory  are  ongoing.  Analytic  second  derivatives  coupled 
with  the  first  derivatives  presented  here  will  provide  the  most 
efficient  route  to  explore  potential  energy  surfaces.  They  also 
provide  a route  to  the  evaluation  of  higher  derivatives. 


APPENDIX:  SIMPLIFICATION  OF  THE  FOURTH-ORDER  FIRST-DERIVATIVE 


Consider  the  fourth-order  derivative  terms  of  the  form  Q(A,BxC); 
(1)  the  quadratic  contribution  to  the  fourth-order  energy  is  given  by 
Q(I,t2(l)xt2(l))/D®,  (2)  the  contributions  to  given  by 
Q(t2(l),IxX2(l))  + Q(t2(l),X2xt2(l>),  and  (3)  the  contribution  to 
r(ij,ab)  of  the  form  Q( t2(l) , t2(l)xt2(l)) • If  one  recalls  that  t2(l) 

= X2(l),  contribution  does  not  need  to  be  computed  for  fourth  order 
perturbation  theory.  This  follows  from  explicitly  considering  the 
relationship  of  each  term  of  the  form  Q(A,BxC). 

The  following  definitions  will  be  used: 


rs 

I = <pq| I rs> 

pq 

ab  ab 

t = t (1)  = <ab||ij>/  (e  + e - e - e ) 
ij  ij  i j a b 

ab 

D =:(c+€  — 6 — c)  • 

ij  i j a b 

Consider  all  contributions  arising  from  the  first  term  of  the 
function  Q(A,BxC)  as  defined  in  Table  (10)  in  the  main  text: 


cd  ab  cd  ab 

I t t / D I 

kl  kl  ij  ij 

cd  ab  cd  ab 

t t I / D II 
kl  kl  ij  ij 
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cd 

ab 

ab  ab 

I 

t / D 

III 

kl 

kl 

kl  ij 

cd 

ab 

cd 

t 

t 

IV 

kl 

kl 

ij 

Now  factorize  term  I minus  term  III  and  term  IV  minus  term  II  and 
compare: 

(I  - III) 


<kl 1 I cd>  <kl 1 I ab>  <kl | | cd>  <kl | | ab>  ) cd  ab 
- t / D 


cd  ab  ' 

<kl| |cd>  <kl| |ab>  t / D 

ij  ij 


ab  cd 

^ ° ^ ° kl  kl 


cd  ab  cd  ab 

' . I < 'a  ♦ 'b  > - ‘ 'c  * '^d  > 1 ' “ , 

kl  kl  ij  ij 


(II-  IV) 


' <ij||cd> 

. (e.+c.-e^-Cd) 


cd  ab 

t t <ij I |cd> 
I kl  kl 


<ij| |cd> 


''  cd  ab 
t t 


(Ei+e.-E^-Eb) 


kl  kl 
cd  ab 


< V ^d^  - / D D_ 


ij  iJ 


cd  ab  cd 


ab 


V,  'u,  ‘ ‘ ^ ‘ ^ “n 

kl  kl  ij  ij 


The  resulting  relation  is  I-III  = IV  - II  or  equivalently,  I + IV  - II 
+ III.  Note  that  the  normal  quadratic  energy  contribution  pairs  with 
the  term  arising  from  Q(t,txt). 
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Next,  consider  the  second  term  of  the  function  Q(A,BxC): 
cd  ac  bd 

-1/2  ABC. 
k.1  ij  kl 

Ignoring  the  factor,  this  term  gives  rise  to  the  following: 


cd 

bd 

ac 

ab 

t 

t 

/D 

I 

kl 

kl  ij 

i ij 

cd 

bd 

ac 

ab 

t 

I 

/D 

II 

kl 

kl 

ij 

ij 

cd 

bd 

ac 

ab 

I 

t 

/D 

III 

kl 

kl 

ij 

ij 

cd 

bd 

ac 

t 

t 

IV 

kl 

kl 

ij 

Now  it  is  necessary  to  factor  terms  I-II  and  IV-III.  The  result 
is  given  explicitly  below: 


(I  - II) 


' <kl| |cd>  <ij I |ac> 

, (e.+e.-e  -e  ) 

' 1 j a c'' 


<kl||cd>  <ij||ac>  bd  ab 

t / D 

^VW^d^  kl  ij 


cd  ac  bd 

t t t [ e.  + e 
kl  ij  kl 

(IV  - III) 

<kl I I bd>  <kl I I bd> 

(Ei+Ej-Ea-Eb) 


ab 

- c_i  - Cj  - ej  . , 


cd  ab 
t t 
kl  kl 
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cd  ac  bd 

kl  ij  kl  ^ iJ 

Here  the  result  is  I-II  + IV-III  = 0 or  again,  I+IV  = II+III.  The 
terms  arising  from  Q(I,txt)  are  matched  with  Q(t,txt)  as  in  the  first 
expression.  In  fact,  the  remaining  terms  factor  in  the  same  manner. 
The  algebraic  manipulation  is  essentially  the  same;  therefore,  only 
the  resulting  expressions  and  relationships  will  be  given  below. 

The  terms  arising  from  the  third  term  of  Q(A,BxC)  are: 


cd  ac  bd  ab 
I t t /D  I 

kl  kl  ij  ij 

cd  ac  bd  ab 

t t I /D  II 

kl  kl  ij  ij 

cd  ac  bd  ab 

t I t /D  III 

kl  kl  ij  ij 

cd  ac  bd 

t t t . IV 

kl  kl  ij 


Considering  I-II  and  III-IV: 


I-II 

cd  bd  ac 
t t t 
kl  ij  kl 

and  III-IV 


ab 

) 

ij 


cd  bd  ac 

' -'k  - "l  * "c  ♦ 'i  * 'j  - "bl  ^ “ij 

kl  13  kl 


The  relationship  I+IV  = II+III  holds  for  this  case,  as  well. 
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Next  consider  term  four  of  the  function  Q(A,BxC): 


cd 

ab 

cd 

ab 

I 

t 

t 

/D 

I 

kl 

ik 

jl 

ij 

cd 

ab 

cd 

ab 

t 

t 

I 

/D 

II 

kl 

ik 

jl 

ij 

cd 

ab 

cd 

ab 

t 

I 

t 

/D 

III 

kl 

ik 

jl 

ij 

cd 

ab 

cd 

t 

t 

t 

• 

IV 

kl 

ik 

jl 

Here,  use  I-III  and  IV-II: 


I-III 


ab  cd  cd 
t t t 
ik  kl  jl 


- 


e 

c 


*=1 


ab 

E,  + e.  ] / D , 
d b 


IV-II 


ab  cd  cd 
t t t [ -E, 

ik  kl  jl 


"l 


+ 


E + 
C 


E.  + 
1 


"d 


/ D 


ab 

ij 


By  this  point,  it  should  be  clear  that  although  no  simple  one-one 
correspondence  exists  between  the  four  Q(A,BxC)  functions,  a secondary 
mapping  occurs  of  the  form  Q(I,txt)  + Q(t,txt)  = Q(t,Ixt)  + Q(t,txl). 
For  completeness,  the  remaining  terms  are  shown  to  obey  this 
relationship. 

Term  five  generates  the  following: 


cd  cd  ab  ab 

I t t /D 
kl  ik  jl  ij 


I 
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cd 

cd 

ab 

ab 

t 

I 

/D 

II 

kl 

ik 

jl 

ij 

cd 

cd 

ab 

ab 

I 

t 

/D 

III 

kl 

ik 

jl 

ij 

cd 

cd 

ab 

t 

t 

IV 

kl 

ik 

jl 

where  one  may  factor  I-II  and  IV-II  to  obtain: 


I-II 


cd  cd 

ab 

ab 

t t 

ik  kl 

t 

jl 

1 \ - 'c  - 'j  - "d  * S'  ^ “jj 

IV-III 

cd  cd 

ab 

ab 

t t 

ik  kl 

t 

jl 

I -s  - S * S * 'j  * s - s'  ' “jj 

Here  I+IV  = 

II+III,  as  expected. 

The 

family  of  expressions  from  term  six, 

cd  ac 

bd 

ab 

I t 

t 

/D  I 

kl  ik 

jl 

ij 

cd  ac 

bd 

ab 

t t 

I 

/D  II 

kl  ik 

jl 

ij 

cd  ac 

bd 

ab 

t I 

t 

/D  III 

kl  ik 

jl 

ij 

cd  ac 

bd 

t t 

t 

IV 

kl  ik 

jl 

provide  similar  results,  if  one  factors  I-II  and  IV-III: 


Ill 


I-II 


ac 

cd 

bd 

ab 

t 

ik 

t 

kl 

t 

jl 

1 "k  - 

e - c.  + e.  1 / D , 

c J b .. 

IV-III 

ac 

cd 

bd 

ab 

t 

t 

t 

+ e + e.  - 1 

C 1 

% 1 ' 

ik 

kl 

jl 

j 

ij 

The 

remaining  term 

follows  the 

pattern: 

cd 

bd 

ac 

ab 

I 

t 

t 

/D 

I 

kl 

ik 

jl 

ij 

cd 

bd 

ac 

ab 

t 

t 

I 

/D 

II 

kl 

ik 

jl 

ij 

cd 

bd 

ac 

ab 

t 

I 

t 

/D 

III 

kl 

ik 

jl 

ij 

cd 

bd 

ac 

t 

t 

t 

« 

IV 

kl 

ik 

jl 

Here,  choose  I-II 

and  IV-III. 

The  result 

summarized  below: 


I-II 


bd 

cd 

ac 

t 

t 

t 

ik 

kl 

jl 

IV-III 

bd 

cd 

ac 

t 

t 

t 

ik 

kl 

jl 

ab 


ab 

D 

ij 


At  this  point,  the  relationship  is  explicit.  Although,  through  the 
course  of  this  exercise,  is  expected  that  the  relationship  will  hold 
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for  each  set  of  terms,  it  is  by  no  means  obvious.  The  importance  of 
this  derivation  is  far-reaching  in  the  development  of  efficient  and 
cost-effective  first  derivatives  of  Q-MBPT(4).  In  addition  to 
reducing  the  number  of  times  the  function  Q(A,BxC)  (which  is  N in  the 
number  of  the  basis  functions)  must  be  computed,  the  resulting 
expressions  follow  the  same  pattern  as  lower  orders  of  MBPT  first 
derivatives.  This  allows  the  computer  codes  that  actually  compute  the 
gradient  expressions  to  function  independently  of  the  level  of 
correlation. 


( 
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